1.1 Le nom des fichiers de codes

Appeler vos fichiers de code avec l’extension .R et donner des noms clairs à vos fichiers de codes, en utilisant éventuellement le symbole underscore ou le trait d’union entre les différents mots. Par exemple :

traitement-base.R
fonctions-utiles.R
exploration.R
prediction.R

1.2 Le nom des objets

Eviter d’appeler vos objets avec des noms de fonctions existants ou de mots clés. Le nombre de caratères devrait être restreint même si le nom doit avoir un sens pour le lecteur. On peut utiliser le symbole underscore ou le point pour séparer des mots entre eux :

arbre_reg
rf_reg
glm_reg

1.3 Espace entre les opérateurs

Aérer au maximum les opérations (affectation, calcul, utilisation de fonctions, etc.) par des espaces :

a <- c(5, NA, 4, 3) 
mean(a, na.rm = TRUE)

Exception faite pour les opérateurs :, :: et ::: où il ne faut pas d’espace :

1:10
stats::lm

1.4 Les conditions

Quand on utilise les conditions if/else, for, while :

• laisser un espace après le mot clé,

• l’accolade ouvrante se trouve à la fin de la ligne contenant le mot clé, l’accolande fermante se trouve sur une nouvelle ligne,

• mettre deux espaces en début de ligne (indentation) à partir de la seconde ligne jusqu’à la fin de la condition,

• si des conditions sont imbriquées, la nouvelle condition devrait se trouver sur la même ligne que l’accolade fermante de la première condition.

if (y < 0 && debug) {
  message("Y is negative")
}

if (y == 0) {
  log(x)
} else {
  y ^ x
}

1.5 La taille d’une ligne

Essayer de limiter 80 caractères à une ligne de code.

1.6 Affectation

Pour affecter une valeur à un objet essayer d’utiliser l’opérateur <- plutôt que =.

a <- 5

Remarque : sur RStudio quelques-une des règles ci-dessus sont implémentées par défaut, notamment les indentations après les boucles ou les fonctions.

Exercice 2.1

Mettre le code suivant en utilisant les règles de style présentées ci-dessus:

my_mean=function(x)
{
# verification
if(!is.numeric(x))
  stop("x must be a numeric vector")
# initialization
n= length(x) 
res =0
for(k in 1:n) 
{
  res= res+x[k]
}
# return res
  return(res/ n)
}  

3.1 Pour mesurer les temps de calcul efficacement

Si vous répétez les instructions précédentes successivement, vous constaterez que les temps de calculs sont à chaque fois différents. Ceci s’explique par le fait que les machines sur lesquelles on travaille exécutent plusieurs tâches à la fois, qui ne sont pas toujours les mêmes à l’instant \(t\). Aussi, une façon de rendre robuste ces temps de calcul est de répéter un certain nombre de fois ces mêmes commandes et de présenter ensuite les résultats statistiques sur ces temps de calcul.

Pour cela, on peut bien évidemment utiliser une boucle for ainsi et faire ensuite des statistiques de base

my_time <- data.frame(time_1 = numeric(10),
                      time_2 = numeric(10))
for (b in 1:10) {
 my_time$time_1[b] <- system.time(trunc_rnorm.1(10000))[3]
 my_time$time_2[b] <- system.time(trunc_rnorm.2(10000))[3]
}
summary(my_time)

##      time_1           time_2       
##  Min.   :0.1100   Min.   :0.02600  
##  1st Qu.:0.1153   1st Qu.:0.02725  
##  Median :0.1185   Median :0.02850  
##  Mean   :0.1185   Mean   :0.02880  
##  3rd Qu.:0.1210   3rd Qu.:0.02975  
##  Max.   :0.1250   Max.   :0.03400

Sinon, on utilise la fonction microbenchmark() du package microbenchmark (Mersmann et al., 2019). En général, on regarde la moyenne.

mbm <- microbenchmark::microbenchmark(
  trunc_rnorm.1(n = 10000),
  trunc_rnorm.2(n = 10000),
  times = 10L)
mbm

## Unit: milliseconds
##                      expr       min        lq      mean    median        uq
##  trunc_rnorm.1(n = 10000) 111.43848 113.39500 118.93426 116.61849 120.07285
##  trunc_rnorm.2(n = 10000)  26.79778  28.71063  30.23288  30.60795  31.59715
##        max neval
##  141.24738    10
##   32.95456    10

Par ailleurs, il est possible de représenter graphiquement ces résultats avec la fonction autoplot() du package ggplot2.

ggplot2::autoplot(mbm)

3.2 Pour comprendre la gestion de la mémoire

Lorsqu’on crée un vecteur sous R, celui-ci est stocké sur un espace mémoire que l’on peut identifier avec la fonction address() du package pryr (Wickham, 2019).

Par exemple :

library("pryr")
x <- numeric(10)
address(x)

## [1] "0x555911a74e78"

Lorsqu’on attribue de nouvelles valeurs à ce vecteur, l’adresse ne change pas :

for (i in 1:10) {
  x[i] <- ifelse(rnorm(1) > 0, 1, 0)
  print(address(x))
}

## [1] "0x55590e6cc628"
## [1] "0x55590e6cc628"
## [1] "0x55590e6cc628"
## [1] "0x55590e6cc628"
## [1] "0x55590e6cc628"
## [1] "0x55590e6cc628"
## [1] "0x55590e6cc628"
## [1] "0x55590e6cc628"
## [1] "0x55590e6cc628"
## [1] "0x55590e6cc628"

En revanche, si on change la taille de ce vecteur, R va systématiquement réserver un nouvel espace mémoire, ce qui est une opération coûteuse en temps calcul, d’où des temps de calcul plus long. Par exemple :

for (i in 11:20) {
  x[i] <- ifelse(rnorm(1) > 0, 1, 0)
  print(address(x))
}

## [1] "0x5559107cba68"
## [1] "0x55591064a078"
## [1] "0x55591064a548"
## [1] "0x55591064a968"
## [1] "0x55591064b3b8"
## [1] "0x55591064b468"
## [1] "0x55590cfdbae0"
## [1] "0x555911c794a0"
## [1] "0x555911c5ce70"
## [1] "0x55590bdc8590"

Exercice 2.2 Comparer les temps de calcul des trois expressions suivantes. Représenter graphiquement la variation des temps de calcul.

n <- 10 ^ 6
# expression 1
x <- numeric(n)
for (k in 1:n) 
  x[k] <- (5 == sample(1:10, 1))
mean(x)
# expression 2
x <- NULL
for (k in 1:n) 
  x <- c(x, (5 == sample(1:10, 1)))
mean(x)
# expression 3
x <- 0
for (k in 1:n)
  x <- x + (5 == sample(1:10, 1))
mean(x)

7.1 apply() pour les matrix/array

Une matrix est une forme particulière d’un array. Une matrix posssède deux dimensions (espace des lignes et espace des colonnes), un array peut avoir autant de dimensions que souhaitées. Les \(D\)-dimensions sont données sous forme d’un vecteur de taille \(D\), chaque élément \(i\) donnant la taille de la dimension \(i\). En général, on utilise peu les array à plus de 3 dimensions. Une caractéristique d’un array est que tous ses éléments sont du même type : numeric, integer, logical, character (rarement complex).

Commençons par créer un array à trois dimensions de taille \(3 \times 2 \times 2\). Pour simplifier la compréhension d’un array, on va donner des noms à chaque composante de chaque dimension. La 1ère dimension correspondra aux nom des individus, la 2ème celle des variables, la 3ème celle des années :

tab <- array(c(100, 90, 110, 25, 24, 28, 175, 180, 190, 68, 74, 85), 
             dim = c(3, 2, 2))
dimnames(tab)[[1]] <- c("Luc", "Pierre", "Pedro")
dimnames(tab)[[2]] <- c("taille", "poids")
dimnames(tab)[[3]] <- c("2000", "2010")
tab

## , , 2000
## 
##        taille poids
## Luc       100    25
## Pierre     90    24
## Pedro     110    28
## 
## , , 2010
## 
##        taille poids
## Luc       175    68
## Pierre    180    74
## Pedro     190    85

La fonction apply() permet de faire des calculs :

• sur une seule dimension. On précise alors le numéro de la dimension dans le deuxième argument et le troisième argument correspond à la fonction qu’on souhaite appliquer sur chaque élément appartenant à cette dimension. Dans ce cas, le résultat est un vecteur de longueur la taille de la dimension choisie. Par exemple, si on souhaite faire la moyenne des éléments en considérant la 2ème dimension.

apply(tab, 2, mean)

##    taille     poids 
## 140.83333  50.66667

Remarque : en réalité, selon la fonction qu’on applique à apply(), le résultat n’est pas forcément un vecteur. Essayer par exemple de remplacer la fonction mean() par range() dans l’exemple ci-dessus.

Equivalent avec les boucles : si on avait utilisé les boucles plutôt que la fonction apply(), on aurait utilisé le code suivant :

res <- numeric(dim(tab)[2])
for (k in seq_along(res))
  res[k] <- mean(tab[, k, ])
res

## [1] 140.83333  50.66667

• sur deux dimensions :

apply(tab, c(2, 3), mean)

##             2000      2010
## taille 100.00000 181.66667
## poids   25.66667  75.66667

Equivalent avec les boucles : si on avait utilisé les boucles pour faire ce calcul, on aurait fait :

res <- array(0, dim = dim(tab)[c(2, 3)])
for (i in 1:dim(tab)[2])
  for (j in 1:dim(tab)[3])
    res[i, j] <- mean(tab[, i, j])
res

##           [,1]      [,2]
## [1,] 100.00000 181.66667
## [2,]  25.66667  75.66667

7.2 Les fonctions colSums(), rowSums(), colMeans(), rowMeans()

Ces fonctions sont équivalentes à la fonction apply() en utilisant FUN=sum ou FUN=mean et en appliquant les bonnes dimensions (1 pour lignes, 2 pour colonnes), mais sont plus rapides car elles ont été codées en langage compilé.

Pour comparer les temps de calcul :

x <- matrix(runif(10e6), nc = 5)
 microbenchmark::microbenchmark(
   apply(x, 2, mean),
   colMeans(x),
   times = 10L)

## Unit: milliseconds
##               expr      min        lq       mean     median         uq
##  apply(x, 2, mean) 98.79311 103.48479 131.653359 111.286857 181.312863
##        colMeans(x)  7.93473   8.31074   8.652441   8.695388   8.888168
##         max neval
##  197.150815    10
##    9.406257    10

7.3 La fonction lapply()

Même si la fonction lapply() s’applique aussi bien sur les vecteurs que les list, on présentera ici son utilisation pour les list/data.frame. L’intérêt de l’appliquer sur des vecteurs sera présenté avec la fonction replicate().

On rappelle ici qu’une list est caractérisée par un certain nombre d’éléments pouvant être de types différents. On accède aux éléments d’une liste par le symbole $ suivi du nom de l’élément de la liste et s’il n’y a pas de nom, on y accède avec les doubles crochets avec l’indice de l’élément auquel on souhaite accèder. Un data.frame est une forme particulière d’une list dans la mesure où les éléments peuvent être vus comme étant les variables (on y accède avec le symbole $).

La fonction lapply() consiste à appliquer la même opération sur chaque élément de la list. Cette opération pouvant être donnée par une fonction de base de R ou bien pouvant être une fonction programmée.

On considère le jeu de données mtcars accessible par défaut sous R et contenant un certaine nombre de variables sur des moteurs de voitures de différentes marques. Pour connaître les caractéristiques moyennes, on applique la fonction mean() à chaque élément de mtcars, autrement dit à chaque variable. Le résultat est retourné sous forme de list de même longeur que la list de départ.

Ici, on a utilisé la fonction unlist() sur le résultat retourné afin de présenter le résultat sous forme d’un vecteur, plus lisible à lire :

unlist(lapply(mtcars, mean))

##        mpg        cyl       disp         hp       drat         wt       qsec 
##  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250  17.848750 
##         vs         am       gear       carb 
##   0.437500   0.406250   3.687500   2.812500

Cette fonction nous permet ainsi d’éviter de faire la boucle suivante :

res <- numeric(length(mtcars))
for (k in seq_along(res))
  res[k] <- mean(mtcars[[k]])
res

##  [1]  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250
##  [7]  17.848750   0.437500   0.406250   3.687500   2.812500

7.4 La fonction sapply()

On a vu que la fonction lapply() retournait une list. Or, si on sait par avance que le résultat retourné pour chaque élément de la liste sera identique pour chaque élément, la fonction sapply() va concaténer les résultats de chaque élément sous forme d’un vecteur ou d’un tableau. Si on reprend l’exemple précédent :

sapply(mtcars, mean)

##        mpg        cyl       disp         hp       drat         wt       qsec 
##  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250  17.848750 
##         vs         am       gear       carb 
##   0.437500   0.406250   3.687500   2.812500

On présente ici un autre exemple où on va créer notre propre fonction qui consiste à retourner, le minimum, le maximum, la moyenne, l’écart-type et la moyenne de chaque élément. On code la fonction suivante où l’argument x sera interpété comme représentant un élément de la list.

f <- function(x) 
  c(min = min(x), max = max(x), mean = mean(x), med = median(x), sd = sd(x))

On applique ensuite la fonction f() à chaque élément de la list avec la fonction sapply(). On transpose le résultat car c’est en général ainsi que sont présentées les statistiques descriptives des variables :

t(sapply(mtcars, f))

##         min     max       mean     med          sd
## mpg  10.400  33.900  20.090625  19.200   6.0269481
## cyl   4.000   8.000   6.187500   6.000   1.7859216
## disp 71.100 472.000 230.721875 196.300 123.9386938
## hp   52.000 335.000 146.687500 123.000  68.5628685
## drat  2.760   4.930   3.596563   3.695   0.5346787
## wt    1.513   5.424   3.217250   3.325   0.9784574
## qsec 14.500  22.900  17.848750  17.710   1.7869432
## vs    0.000   1.000   0.437500   0.000   0.5040161
## am    0.000   1.000   0.406250   0.000   0.4989909
## gear  3.000   5.000   3.687500   4.000   0.7378041
## carb  1.000   8.000   2.812500   2.000   1.6152000

Par ailleurs, pour insérer une table dans un document Markdown, on peut utiliser la fonction kbl() du package kableExtra dont on verra d’autres exemples dans le chapitre sur la visualisation de données.

kableExtra::kbl(t(round(sapply(mtcars, f), 3))) 

Remarque : si on on souhaite appliquer la fonction f() pour les éléments qui sont de type numeric et la fonction table() pour les éléments de type factor, alors la fonction sapply() retournera le même résultat que la fonction lapply(), à savoir une list, car les résultats retournés pour chaque élément ne seront pas les mêmes :

sapply(iris, function(x){
  if (is.numeric(x)) 
    f(x)
  else
    table(x)
})

## $Sepal.Length
##       min       max      mean       med        sd 
## 4.3000000 7.9000000 5.8433333 5.8000000 0.8280661 
## 
## $Sepal.Width
##       min       max      mean       med        sd 
## 2.0000000 4.4000000 3.0573333 3.0000000 0.4358663 
## 
## $Petal.Length
##      min      max     mean      med       sd 
## 1.000000 6.900000 3.758000 4.350000 1.765298 
## 
## $Petal.Width
##       min       max      mean       med        sd 
## 0.1000000 2.5000000 1.1993333 1.3000000 0.7622377 
## 
## $Species
## x
##     setosa versicolor  virginica 
##         50         50         50

7.5 La fonction replicate()

On va partir de l’exemple suivant : on souhaitre simuler 5 échantillons de taille 10 distribués selon une loi uniforme \(U_{[0,1]}\).

• La première solution serait de faire une boucle :

res <- vector("list", 5)
for (k in 1:5)
  res[[k]] <- runif(10)

• La deuxième solution consiste à appliquer la fonction lapply() sur un vecteur quelconque de taille 5 (ici, le plus simple est de faire 1:5, mais on aurait pu prendre n’importe quel vecteur de taille 5) et de lui appliquer une fonction qui retourne un vecteur simulé selon une \(U_{[0,1]}\). Un vecteur est donc considéré ici comme une liste où chaque élément du vecteur serait un élement d’une liste.

res <- lapply(1:5, function(x) runif(10))

Remarque : dans le deuxième argument, la fonction qu’on va appliquer sur chaque élément d’une vecteur prend comme argument d’entrée x (un élément du vecteur), mais on ne l’utilise pas à l’intérieur de la fonction, car on n’en a pas besoin pour faire ce qu’on souhaite faire.

• La troisième solution consiste à utiliser la fonction replicate(), qui comme la fonction sapply() va retourner le résultat sous une forme plus simplifiée qu’une list :

res <- replicate(5, runif(10))

Remarque : la fonction replicate() fait précisément appel à la fonction sapply(). Pour cela, elle créé un vecteur de 0L de taille le premier argument de replicate().

7.6 La fonction mapply()

Pour illustrer cette fonction, on considère 5 échantillons de taille 10 issus d’une loi \(U_{[0,1]}\). Pour calculer la moyenne par échantilon, on peut utiliser la fonction lapply() :

xs <- replicate(5, runif(10), simplify = FALSE)
lapply(xs, mean)

## [[1]]
## [1] 0.4905019
## 
## [[2]]
## [1] 0.493022
## 
## [[3]]
## [1] 0.4635896
## 
## [[4]]
## [1] 0.4134804
## 
## [[5]]
## [1] 0.3997965

A présent, supposons qu’on souhaite pondérer les moyennes par un vecteur de poids qui est différent selon chaque échantillon. On créé ici les poids associés à chaque échantillon sous forme d’une liste :

ws <- replicate(5, rpois(10, 5) + 1, simplify = FALSE)

Une façon de faire avec des boucles serait la suivante :

res <- numeric(length(xs))
for (k in seq_along(res))
  res[k] <- sum(xs[[k]] * ws[[k]]) / sum(ws[[k]])
res

## [1] 0.5105359 0.4561600 0.5012219 0.4297210 0.3654479

Le problème de la fonction lapply() est qu’on ne peut pas appliquer une fonction sur deux list simultanément. C’est ce que fait la fonction mapply() qui prend comme 1er argument d’entrée une fonction à appliquer sur les éléments d’autant de list que l’on souhaite. La fonction qu’on applique contient donc deux arguments d’entrée x et y qui correspondent aux éléments de chacune des deux listes :

mapply(function(x, y) sum(x * y) / mean(y), xs, ws)

## [1] 5.105359 4.561600 5.012219 4.297210 3.654479

7.7 La fonction tapply()

On considère les deux vecteurs suivants :

(taille <- c(rnorm(5, 165, 10), rnorm(5, 175, 10)))

##  [1] 177.2044 181.6881 167.4184 179.6390 170.8696 161.6642 184.5007 181.4993
##  [9] 171.5086 172.9919

(sexe <- rep(c("M", "F"), each = 5))

##  [1] "M" "M" "M" "M" "M" "F" "F" "F" "F" "F"

On souhaite calculer la moyenne dans chaque sous-groupe “M” et “F”. Une façon de faire est de créer une list contenant deux éléments, le premier correspondant au vecteur de taille du sous-groupe “M” et le second correspondant au vecteur de taille du sous-groupe “F”. C’est ce que fait la fonction split() ci-dessous :

(res.split <- split(taille, sexe))

## $F
## [1] 161.6642 184.5007 181.4993 171.5086 172.9919
## 
## $M
## [1] 177.2044 181.6881 167.4184 179.6390 170.8696

Ensuite, il suffit d’appliquer la fonction sapply() à cette list :

sapply(res.split, mean)

##        F        M 
## 174.4329 175.3639

La fonction tapply() permet de faire ce calcul en une seule ligne :

tapply(taille, sexe, mean)

##        F        M 
## 174.4329 175.3639

Remarque : la fonction tapply() ne s’appliquant que sur des vecteurs, la fonction by() permet de généraliser la fonction tapply() aux data.frame.

Exercise 2.5

• Créer deux fonctions mean_rnorm.1() et mean_rnorm.2() qui prennent en entrée les paramètres n et p. Ces fonctions permettent de simuler n échantillons de taille p de lois normales centrées réduites et retournent pour chaque échantillon la moyenne. mean_rnorm.1() utilisera des boucles, mean_rnorm.2() la fonction replicate().

• Comparer les temps de calcul de ces 2 fonctions.

8.1 Créer des sous-fonctions

Il ne faut pas hésiter à créer des petites fonctions qui peuvent être

• locales si elles ne sont utilisées qu’à l’intérieur d’une seule fonction principale

• globales si elles ont à vocation d’être uilisées dans plusieurs fonctions.

Supposons par exemple qu’on souhaite calculer une densité non paramétrique de la densité en utilisant un des noyaux suivants :

• noyau biweight \(K(x) = \frac{15}{16}(1-(\frac{x}{h})^2)^21_{(\frac{x}{h})^2\leq 1}\)

• noyau triweight \(K(x) = \frac{35}{32}(1-(\frac{x}{h})^2)^31_{(\frac{x}{h})^2\leq 1}\)

• noyau gaussien \(K(x) = \frac{1}{\sqrt{2\pi}}\exp(-0.5(\frac{x}{h})^2))\)

où h est le paramètre de lissage. Dans un premier temps, pn peut alors programmer une fonction par noyau :

biweight <- function(x, h) {
  return(15 / 16 * (1 - (x / h) ^ 2) ^ 2 * ifelse((x / h) ^ 2 <= 1, 1, 0))
}
triweight <- function(x, h) {
  return(35 / 32*(1 - (x / h) ^ 2) ^ 3 * ifelse((x / h) ^ 2 <= 1, 1, 0))
}
gaussian <- function(x, h) {
  return(1 / sqrt(2 * pi) * exp(-0.5 * (x / h) ^ 2))
}

8.2 Structure de contrôle

Les structures de contrôle classiques (for, while, repeat, if/else) sont bien entendu disponibles dans R (voir les document “Pour se donner un peu d’R” et “Introduction à R”. Nous nous intéressons ici à la notion d’aiguillage (fonction switch()) moins couramment utilisée bien que très pratique.

Celle-ci permet d’éviter d’emboîter des if/else lorsqu’on a plusieurs options possibles. Supposons par exemple qu’on souhaite calculer une densité non paramétrique de la densité en utilisant un des noyaux suivants :

• noyau biweight \(K(x) = \frac{15}{16}(1-(\frac{x}{h})^2)^21_{(\frac{x}{h})^2\leq 1}\)

• noyau triweight \(K(x) = \frac{35}{32}(1-(\frac{x}{h})^2)^31_{(\frac{x}{h})^2\leq 1}\)

• noyau gaussien \(K(x) = \frac{1}{\sqrt{2\pi}}\exp(-0.5(\frac{x}{h})^2))\)

où h est le paramètre de lissage. On peut alors programmer sous R la fonction f_noyau() suivante, prenant en entrée les arguments x, h et type où type est un caractère prenant les valeurs “bi”, “tri” et “gauss”.

f_noyau <- function(x, h, type = "bi") {
  if (type == "bi") {
    biweight(x, h)
  } else {
    if (type == "tri") {
       triweight(x, h)
    } else {
      gaussian(x, h)
    }
  }
}

f_noyau.2 <- function(x, h, type = "bi") {
  switch(type, bi = biweight(x, h),
               tri = triweight(x, h),
               gauss = gaussian(x, h),
               "Préciser une autre valeur de type")
}

x <- seq(-1, 1, 0.01)
plot(x, f_noyau.2(x, 0.3, type = "bi"), type = "l", ylab = "", ylim = c(0, 1.2))
lines(x, f_noyau.2(x, 0.3, type = "tri"), lty = 2)
lines(x, f_noyau.2(x, 0.3, type = "gauss"), lty = 3)

f_noyau.2(x, 0.3, type = "unif")

## [1] "Préciser une autre valeur de type"

8.3 Fonction stopifnot()

On connaît la fonction stop() qui s’insère en général à l’intérieur des conditions if/else.

La fonction stopifnot() permet de tester simultanément plusieurs conditions et évite donc d’avoir recours aux if/else. Elle est particulièrement utile à l’intérieur d’une fonction pour, par exemple, vérifier la conformité des valeurs passées en paramètres. Un exemple d’uitlisation :

stopifnot(1 < 2, length(1:2) == 2, pi < 2, cos(pi) > 3)

## Error: pi < 2 n'est pas TRUE

Si plusieurs conditions ne sont pas respectées, c’est la première non respectée qui est renvoyée comme source de l’erreur.

Exercise 2.6

Modifier la fonction f_noyau.2() définie précédemment pour vérifier que :

• le paramètre x est numeric

• h est de taille 1 et de type numeric

• type est une chaîne de caractères (on supprimera la possibilité de passer un entier pour choisir le type).

8.4 Les arguments dans une fonction

f_1 <- function(a = 5, b = 4, d = 3, e = 1)
  (a + b)^2

f_2 <- function(a = 5, b = 4)
  (a + b)^2

microbenchmark::microbenchmark(
  f_1(),
  f_2()
)

## Unit: nanoseconds
##   expr min    lq     mean median  uq     max neval
##  f_1() 307 318.5 12118.81  337.5 356 1173445   100
##  f_2() 243 253.0 11388.95  269.5 285 1109342   100

randomise <- function(f) f(runif(1e3))

randomise(mean)

## [1] 0.4929316

randomise(mean)

## [1] 0.4926691

randomise(sum)

## [1] 486.2975

f_power <- function(exponent)
  function(x) x^exponent 

f_power(2)(1:5)

## [1]  1  4  9 16 25

f_power(3)(1:5)

## [1]   1   8  27  64 125

f_square <- f_power(2)
f_cube <- f_power(3)
class(f_square)

## [1] "function"

class(f_cube)

## [1] "function"

f_square(1:5)

## [1]  1  4  9 16 25

f_cube(1:5)

## [1]   1   8  27  64 125

plot_reg <- function(x, y, np = TRUE, ...) {
  plot(y ~ x, ...)
  abline(lm(y ~ x), col = "blue")
  if (np) {
    np.reg <- loess(y ~ x)
    x.seq <- seq(min(x), max(x), length.out = 25)
    lines(x.seq, predict(np.reg, x.seq), col = "red")
  }
}

x <- rnorm(100) 
y <- x*runif(100)

plot_reg(x, y, np = FALSE)

plot_reg(x, y, pch = 16, col = "pink",
 xlab = "variable explicative", ylab = "variable à expliquer")

9.1 Les fonctions quote() et substitute()

Les fonctions quote() et substitute() permettent de capturer une instruction sans l’évaluer. Pour généraliser, on va appeler cette capture une expression. La différence entre eval() et substitute() est décrite dans le chapitre Non-standard evaluation de l’ouvrage d’Hadley Wickham, mais pour résumer quote() ne fait que retourner l’expression qu’on lui donne en entrée alors que substitute est un peu plus complexe.

quote(1:10)

## 1:10

quote(f(1:10))

## f(1:10)

quote(f(x))

## f(x)

quote(f(x + y^2))

## f(x + y^2)

substitute(1:10)

## 1:10

substitute(f(1:10))

## f(1:10)

substitute(f(x))

## f(x)

substitute(f(x + y^2))

## f(x + y^2)

Une fois qu’on a fait cette capture de code, le but est d’évaluer cette expression dans un environnement spécifique et au moment où on le souhaitera avec la fonction eval().

9.2 La fonction eval()

Cette fonction permet d’évaluer une expression. Par défaut, elle est évaluée dans l’environnement courant de R. Autrement dit, si on considère l’expression suivante :

a <- quote(f.eval(x.eval))

Si on veut évaluer cette expression, il faut que la fonction f.eval() et l’objet x.eval existent tous les deux dans l’environnement où on souhaite exécuter cette expression. Par défaut, eval() ira chercher dans l’environnement courant. Ici, on obtient le message d’erreur suivant car la fonction f.eval() et l’objet x.eval n’ont pas été définies au préalables :

eval(a)

## Error in f.eval(x.eval): impossible de trouver la fonction "f.eval"

On considère un autre exemple avec l’expression suivante :

a <- quote(Species == "setosa")

La commande suivante donnera un message d’erreur :

eval(a)

## Error in eval(a): objet 'Species' introuvable

En revanche, on peut préciser qu’on va trouver le nom Species dans l’environnement du data.frame nommé iris :

eval(a, envir = iris)

##   [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [13]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [25]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [37]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [49]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] FALSE FALSE FALSE FALSE FALSE FALSE

subset(iris, Species == "setosa" & Sepal.Length > 5.5)

##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 15          5.8         4.0          1.2         0.2  setosa
## 16          5.7         4.4          1.5         0.4  setosa
## 19          5.7         3.8          1.7         0.3  setosa

iris[iris$Species == "setosa" & iris$Sepal.Length > 5.5, ]

##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 15          5.8         4.0          1.2         0.2  setosa
## 16          5.7         4.4          1.5         0.4  setosa
## 19          5.7         3.8          1.7         0.3  setosa

plot(1:10, rnorm(10))

9.3 Débugger une fonction

sum(c("a", "b"))

## Error in sum(c("a", "b")): 'type' (character) de l'argument incorrect

ex_bug <- function(n, p) {
  n2 <- as.character(n)
  length(which(mean_rnorm.1(n2, p) > mean_rnorm.2(n, p)))
}

ex_bug(5, 10)

## Error in mean_rnorm.1(n2, p): impossible de trouver la fonction "mean_rnorm.1"

traceback()

4: matrix(0, n, 1) at #2
3: mean_rnorm.1(n2, p) at #3
2: which(mean_rnorm.1(n2, p) > mean_rnorm.2(n, p)) at #3
1: ex_bug(5, 10)

ex_bug.2 <- function(x){
  x <- x + 1
  x <- x^3
  x <- log(x)
x
}

debugonce(ex_bug.2)
ex_bug.2(-5)

f_error.1 <- function(x) {
  try(x <- log(x))
  x
}

f_error.1("10")

## Error in log(x) : argument non numérique pour une fonction mathématique

## [1] "10"

10.1 La norme S3

Nous nous sommes inspirés du document suivant pour décrire la programmation orientée qui utilise la norme S3 : http://www.duclert.org/r-divers/classes-S3-R.php. Pour expliquer son principe, nous allons prende un exemple concret.

On voudrait calculer l’aire de différentes figures géométriques :

• un carré, dont il nous faut connaître la longueur d’un côté, par exemple :

squ <- 3

• un rectangle, défini par sa longueur et sa largeur, par exemple :

rec <- c(5, 6)

• un cercle, caractérisé par son rayon, par exemple :

cir <- sqrt(10)

Pour ce faire, on souhaiterait utiliser une fonction qui s’appelle getArea() et qui puisse reconaître quand on l’applique à un carré, un rectangle ou un cercle afin d’utiliser la bonne formule. Ceci consiste à définir une méthode générique getArea qui renvoie à la bonne implémentation en fonction de la classe de l’argument d’entrée. C’est ce qu’on appelle le polymorphisme.

class(cir) <- "cercle"
cir

## [1] 3.162278
## attr(,"class")
## [1] "cercle"

class(rec) <- "rectangle"
rec

## [1] 5 6
## attr(,"class")
## [1] "rectangle"

class(squ) <- "carre" 
squ

## [1] 3
## attr(,"class")
## [1] "carre"

area <- function(x) {
   switch(class(x), 
      carre = x ^ 2,
      rectangle = x[1] * x[2],
      cercle = pi * x ^ 2,
      "class should be among carre/rec/cercle")
}

getArea <- function(obj)
  UseMethod("getArea", obj)

getArea.default <- function(obj) {
  stop("Méthode getArea non définie pour ce type d'objet")
}

getArea.cercle <- function(obj) {
  pi * obj[1] ^ 2
}

getArea.rectangle <- function(obj) {
  obj[1] * obj[2]
}

getArea.carre <- function(obj) {
  obj[1]^2
}

getArea(cir)

## [1] 31.41593

getArea(rec)

## [1] 30

getArea(squ)

## [1] 9

a <- 5
getArea(a)

## Error in getArea.default(a): Méthode getArea non définie pour ce type d'objet

plot(rnorm(10), runif(10))
plot(lm(Sepal.Length ~ Sepal.Width, data = iris))

library("pryr")

df <- data.frame(x = 1:10, y = letters[1:10])
pryr::otype(df)    

## [1] "S3"

methods("plot")

##  [1] plot,ANY-method     plot,color-method   plot.acf*          
##  [4] plot.data.frame*    plot.decomposed.ts* plot.default       
##  [7] plot.dendrogram*    plot.density*       plot.ecdf          
## [10] plot.factor*        plot.formula*       plot.function      
## [13] plot.ggplot*        plot.gtable*        plot.hcl_palettes* 
## [16] plot.hclust*        plot.histogram*     plot.HoltWinters*  
## [19] plot.isoreg*        plot.lm*            plot.medpolish*    
## [22] plot.mlm*           plot.numpy.ndarray* plot.ppr*          
## [25] plot.prcomp*        plot.princomp*      plot.profile.nls*  
## [28] plot.R6*            plot.raster*        plot.shingle*      
## [31] plot.spec*          plot.stepfun        plot.stl*          
## [34] plot.table*         plot.trans*         plot.trellis*      
## [37] plot.ts             plot.tskernel*      plot.TukeyHSD*     
## see '?methods' for accessing help and source code

methods(class = "lm")

##  [1] add1           alias          anova          case.names     coerce        
##  [6] confint        cooks.distance deviance       dfbeta         dfbetas       
## [11] drop1          dummy.coef     effects        extractAIC     family        
## [16] formula        fortify        hatvalues      influence      initialize    
## [21] kappa          labels         logLik         model.frame    model.matrix  
## [26] nobs           plot           predict        print          proj          
## [31] qr             residuals      rstandard      rstudent       show          
## [36] simulate       slotsFromS3    summary        variable.names vcov          
## see '?methods' for accessing help and source code

10.2 La norme S4

Ici, on présentera uniquement les grandes lignes de la norme S4. Cette section est inspirée du chapitre OO field guide du livre d’Hadley Wickham.

La norme S4 s’inspire de la norme S3 auxquelles s’ajoutent certaines caractéristiques telles que :

• les classes sont définies en ajoutant un certain nombre de règles et de précisions les concernant. Pour mieux comprendre, dans l’exemple ci-dessous, on peut définir l’objet a comme appartenant à la classe lm dans la norme S3, alors qu’il n’en possède pas les caractéristiques :

a <- 5
class(a) <- "lm"

• on utilisera un opérateur spécial @ pour extraire des éléments d’un objet ayant la norme S4.

Pour savoir si un objet appartient à la norme S4, on peut utiliser la fonction isS4(). La plupart du temps, la norme S4 n’a pas été utilisée pour définir des objets et fonctions qui appartiennant à l’environnement de base. Ils sont donc plutôt présents dans des nouveaux packages qu’on va charger.

Les principales fonctions utilisées pour travailler sur la norme S4 sont :

• la fonction setClass() qui permet de définir une nouvelle classe. Le 1er argument est le nom de la classe, le second est le résultat de la fonction representation() qui définit le nom et le type des éléments que la classe contient et le troisième argument est le résultat de la fonction prototype() qui donne des valeurs par défaut. Par exemple :

setClass("Personne", representation(nom = "character", age = "numeric"),
                        prototype(nom = NA_character_, age = NA_real_))

• la fonction new() permet de créer un nouvel objet.

hadley <- new("Personne", nom = "Hadley", age = 31)

Remarque : on accède aux éléments d’un objet de norme S4 avec la commande spéciale @ :

hadley@age

## [1] 31

Les autres fonctions utiles dont vous trouverez des exemples d’utilisation dans http://adv-r.had.co.nz/S4.html sont :

• la fonction setMethods() permettant de définir des nouvelles méthodes.

• les fonctions as() et setAs() pour passer d’une classe d’objet à une autre (quand cela est possible).

• les fonctions setValidity() et validObject() pour vérifier la validité.

• les fonctions showClass(), showMethods() et getMethod() pour accéder aux propriétés des objets et méthodes créées.

Bilbiographie : la présentation de F. Leisch à useR! 2004 (http://www.ci.tuwien.ac.at/Conferences/useR-2004/Keynotes/Leisch.pdf) ainsi que le manuel de C. Genolini (https://cran.r-project.org/doc/contrib/Genolini-PetitManuelDeS4.pdf).

11.1 La fonction est dans l’environnement courant

Vous tapez le nom de la fonction dans la console et le code apparaît :

sapply

## function (X, FUN, ..., simplify = TRUE, USE.NAMES = TRUE) 
## {
##     FUN <- match.fun(FUN)
##     answer <- lapply(X = X, FUN = FUN, ...)
##     if (USE.NAMES && is.character(X) && is.null(names(answer))) 
##         names(answer) <- X
##     if (!isFALSE(simplify)) 
##         simplify2array(answer, higher = (simplify == "array"))
##     else answer
## }
## <bytecode: 0x55590c947910>
## <environment: namespace:base>

11.2 La fonction est une méthode générique de type S3

Vous tapez le nom d’une fonction lambda() et vous obtenez une ligne UseMethod(“lambda”). Ceci signifie que lambda est une méthode générique pour des objets issus de la classe S3. Il existe donc plusieurs fonctions associées à la fonction lambda(). Pour connaître quelles sont les différentes fonctions associées à la fonction lambda(), on utilise la fonction methods(). Par exemple, pour la fonction summary(), si on tape le nom de la fonction dans la console, on obtient :

summary

## function (object, ...) 
## UseMethod("summary")
## <bytecode: 0x55590df4b3a8>
## <environment: namespace:base>

Pour obtenir les différentes fonctions associées à summary(), on utilise donc la fonction methods().

methods("summary")

##  [1] summary,ANY-method                   summary,DBIObject-method            
##  [3] summary,diagonalMatrix-method        summary,sparseMatrix-method         
##  [5] summary.aov                          summary.aovlist*                    
##  [7] summary.aspell*                      summary.check_packages_in_dir*      
##  [9] summary.connection                   summary.data.frame                  
## [11] summary.Date                         summary.default                     
## [13] summary.ecdf*                        summary.factor                      
## [15] summary.ggplot*                      summary.glm                         
## [17] summary.hcl_palettes*                summary.infl*                       
## [19] summary.lm                           summary.loess*                      
## [21] summary.manova                       summary.matrix                      
## [23] summary.microbenchmark*              summary.mlm*                        
## [25] summary.nls*                         summary.packageStatus*              
## [27] summary.pandas.core.frame.DataFrame* summary.pandas.core.series.Series*  
## [29] summary.POSIXct                      summary.POSIXlt                     
## [31] summary.ppr*                         summary.prcomp*                     
## [33] summary.princomp*                    summary.proc_time                   
## [35] summary.python.builtin.object*       summary.rlang_error*                
## [37] summary.rlang_trace*                 summary.shingle*                    
## [39] summary.srcfile                      summary.srcref                      
## [41] summary.stepfun                      summary.stl*                        
## [43] summary.table                        summary.trellis*                    
## [45] summary.tukeysmooth*                 summary.vctrs_sclr*                 
## [47] summary.vctrs_vctr*                  summary.warnings                    
## see '?methods' for accessing help and source code

Pour lire le code des fonctions affichées ci-dessus, il y a deux options :

• si il n’y pas d’étoile à côté de la fonction, vous pouvez obtenir directement le code de la fonction en tapant son nom complet :

summary.proc_time

## function (object, ...) 
## {
##     if (!is.na(object[4L])) 
##         object[1L] <- object[1L] + object[4L]
##     if (!is.na(object[5L])) 
##         object[2L] <- object[2L] + object[5L]
##     object <- object[1L:3L]
##     names(object) <- c(gettext("user"), gettext("system"), gettext("elapsed"))
##     object
## }
## <bytecode: 0x55590df4b808>
## <environment: namespace:base>

• s’il y a une étoile, c’est que la fonction se trouve dans un package et qu’on doit indiquer le nom du package suivi de ::: :

stats:::summary.ecdf

Si on ne connaît pas le nom du package, on peut alors utiliser la fonction getAnywhere() qui va chercher dans toutes les librairies chargées :

getAnywhere(summary.ecdf)

## A single object matching 'summary.ecdf' was found
## It was found in the following places
##   registered S3 method for summary from namespace stats
##   namespace:stats
## with value
## 
## function (object, ...) 
## {
##     n <- length(eval(expression(x), envir = environment(object)))
##     header <- paste("Empirical CDF:\t ", n, "unique values with summary\n")
##     structure(summary(knots(object), ...), header = header, class = "summary.ecdf")
## }
## <bytecode: 0x5559119d7b88>
## <environment: namespace:stats>

11.3 La fonction fait appel à du code C

Si on affiche le code la fonction sum(), on constate que le calcul n’est pas directement fait dans le corps de la fonction. C’est ce que nous avons vu précédemment, à savoir qu’il existe un grand nombre de fonctions qui font appel à du code compilé pour gagner en performances. Lorsqu’on voit dans une fonction les instructions .Primitive() ou .Internal(), cela signifie que R va en réalité exécuter des fonctions qui ont été compilées en C. La différence entre .Primitive() ou .Internal() vient du fait que dans le deuxième cas, on peut appeler des fonctions qui sont codées en langage R.

sum

## function (..., na.rm = FALSE)  .Primitive("sum")

On peut donc faire appel à du code C de plusieurs façons. On verra ici les deux principales façons de le faire.

pryr::show_c_source(.Internal(mean(x)))

qnorm

## function (p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) 
## .Call(C_qnorm, p, mean, sd, lower.tail, log.p)
## <bytecode: 0x555914856098>
## <environment: namespace:stats>

site:https://svn.r-project.org/R/trunk/src qnorm

11.4 La fonction est une méthode générique de type S4

C’est le cas des fonctions programmées dans le package Matrix.

require("Matrix")

Si on utiliser la fonction methods() sur la méthode dim, pour identifier les méthodes génériques qui s’appliquent à des objets S4, on constate que le nom de la fonction est suivi d’une virgule, suivi alors de la classe de l’objet sur laquelle s’applique la méthode. Dans le cas de dim(), elle peut s’appliquer sur les classes d’objet MatrixFactorization ou Matrix.

methods("dim")

## [1] dim,Matrix-method                dim,MatrixFactorization-method  
## [3] dim.data.frame                   dim.gtable*                     
## [5] dim.layout*                      dim.pandas.core.frame.DataFrame*
## [7] dim.pandas.core.series.Series*   dim.scipy.sparse.base.spmatrix* 
## [9] dim.trellis*                    
## see '?methods' for accessing help and source code

Pour obtenir le code source de la fonction dim() qui s’applique à la classe d’objet MatrixFactorization, on utilise la fonction getMethod() ainsi :

getMethod("dim", "MatrixFactorization")

## Method Definition:
## 
## function (x) 
## x@Dim
## <bytecode: 0x55590f4ed988>
## <environment: namespace:Matrix>
## 
## Signatures:
##         x                    
## target  "MatrixFactorization"
## defined "MatrixFactorization"

Pour connaître les méthodes S4 qui s’appliquent sur une classe d’objets particuliers, on utilise la fonction showMethods() :

showMethods(class = "MatrixFactorization")

## 
## Function "asJSON":
##  <not an S4 generic function>
## 
## Function "complete":
##  <not an S4 generic function>
## 
## Function "coords":
##  <not an S4 generic function>
## 
## Function "dbAppendTable":
##  <not an S4 generic function>
## 
## Function "dbBegin":
##  <not an S4 generic function>
## 
## Function "dbBind":
##  <not an S4 generic function>
## 
## Function "dbCallProc":
##  <not an S4 generic function>
## 
## Function "dbCanConnect":
##  <not an S4 generic function>
## 
## Function "dbClearResult":
##  <not an S4 generic function>
## 
## Function "dbColumnInfo":
##  <not an S4 generic function>
## 
## Function "dbCommit":
##  <not an S4 generic function>
## 
## Function "dbConnect":
##  <not an S4 generic function>
## 
## Function "dbCreateTable":
##  <not an S4 generic function>
## 
## Function "dbDataType":
##  <not an S4 generic function>
## 
## Function "dbDisconnect":
##  <not an S4 generic function>
## 
## Function "dbDriver":
##  <not an S4 generic function>
## 
## Function "dbExecute":
##  <not an S4 generic function>
## 
## Function "dbExistsTable":
##  <not an S4 generic function>
## 
## Function "dbFetch":
##  <not an S4 generic function>
## 
## Function "dbGetConnectArgs":
##  <not an S4 generic function>
## 
## Function "dbGetException":
##  <not an S4 generic function>
## 
## Function "dbGetInfo":
##  <not an S4 generic function>
## 
## Function "dbGetQuery":
##  <not an S4 generic function>
## 
## Function "dbGetRowCount":
##  <not an S4 generic function>
## 
## Function "dbGetRowsAffected":
##  <not an S4 generic function>
## 
## Function "dbGetStatement":
##  <not an S4 generic function>
## 
## Function "dbHasCompleted":
##  <not an S4 generic function>
## 
## Function "dbiDataType":
##  <not an S4 generic function>
## 
## Function "dbIsReadOnly":
##  <not an S4 generic function>
## 
## Function "dbIsValid":
##  <not an S4 generic function>
## 
## Function "dbListConnections":
##  <not an S4 generic function>
## 
## Function "dbListFields":
##  <not an S4 generic function>
## 
## Function "dbListObjects":
##  <not an S4 generic function>
## 
## Function "dbListResults":
##  <not an S4 generic function>
## 
## Function "dbListTables":
##  <not an S4 generic function>
## 
## Function "dbQuoteIdentifier":
##  <not an S4 generic function>
## 
## Function "dbQuoteLiteral":
##  <not an S4 generic function>
## 
## Function "dbQuoteString":
##  <not an S4 generic function>
## 
## Function "dbReadTable":
##  <not an S4 generic function>
## 
## Function "dbRemoveTable":
##  <not an S4 generic function>
## 
## Function "dbRollback":
##  <not an S4 generic function>
## 
## Function "dbSendQuery":
##  <not an S4 generic function>
## 
## Function "dbSendStatement":
##  <not an S4 generic function>
## 
## Function "dbSetDataMappings":
##  <not an S4 generic function>
## 
## Function "dbUnloadDriver":
##  <not an S4 generic function>
## 
## Function "dbUnquoteIdentifier":
##  <not an S4 generic function>
## 
## Function "dbWithTransaction":
##  <not an S4 generic function>
## 
## Function "dbWriteTable":
##  <not an S4 generic function>
## Function: dim (package base)
## x="MatrixFactorization"
## 
## Function: expand (package Matrix)
## x="MatrixFactorization"
## 
## 
## Function "fetch":
##  <not an S4 generic function>
## 
## Function "functions":
##  <not an S4 generic function>
## 
## Function "isSQLKeyword":
##  <not an S4 generic function>
## 
## Function "make.db.names":
##  <not an S4 generic function>
## 
## Function "plot":
##  <not an S4 generic function>
## Function: show (package methods)
## object="MatrixFactorization"
## 
## Function: solve (package base)
## a="MatrixFactorization", b="ANY"
## a="MatrixFactorization", b="missing"
## a="MatrixFactorization", b="numeric"
## 
## 
## Function "sqlAppendTable":
##  <not an S4 generic function>
## 
## Function "sqlCreateTable":
##  <not an S4 generic function>
## 
## Function "sqlData":
##  <not an S4 generic function>
## 
## Function "sqlInterpolate":
##  <not an S4 generic function>
## 
## Function "SQLKeywords":
##  <not an S4 generic function>
## 
## Function "sqlParseVariables":
##  <not an S4 generic function>