Formation R perfectionnement Martin Chevalier (Insee)

---

Formation R Perfectionnement

Cette page comporte l’ensemble des cas pratiques de la formation R perfectionnement donnée à la Drees les 16 et 17 avril 2018 accompagnés de leur correction.

• Support de présentation
• Données utilisées dans les cas pratiques

Les supports de cette formation ont été conçus sous RStudio avec R Markdown et compilés le 10/04/2018. Certains éléments de mise en forme du site compagnon sont repris de l’ouvrage R packages de Hadley Wickham.

Ces supports seront durablement disponibles à l’adresse http:// t.slmc.fr/perf et leur code source sur github. L’ensemble est librement réutilisable sous © 2016-2018 Martin Chevalier CC BY-NC-SA 3.0.

Savoir utiliser les fonctions *apply(), do.call() et Reduce()

Le package microbenchmark est souvent utilisé dans cette partie et les suivantes pour mesurer la performance des des solutions testées :

install.packages("microbenchmark")
library(microbenchmark)

 

Cas pratique 1 Fonctions et environnements

Tout ce qui se passe dans R correspond à un appel de fonction. Comprendre le fonctionnement des fonctions et savoir en créer soi-même est donc crucial.

1. Utilisez les guillemets simples inversés (AltGr + 7 sur le clavier azerty) pour afficher le code associé au signe +. Utilisez-le comme une fonction classique avec la syntaxe nomFonction(parametre1, parametre2).

   Afficher/masquer la solution

   ---

   # Pour afficher le code d'une fonction, il suffit
   # de saisir son nom. Dans le cas des signes 
   # arithmétiques, il convient d'entourer le nom
   # de guillemets inversés
   `+`
   ## function (e1, e2)  .Primitive("+")
   
   # En utilisant les guillemets inversés, on peut
   # utiliser la fonction `+`() comme n'importe
   # quelle fonction (en appelant ses arguments dans
   # la parenthèse).
   `+`(2, 2)
   ## [1] 4

   ---

    

2. Définissez la fonction monCalcul(x, puissance) qui pour un vecteur numérique x quelconque :

   1. calcule sa somme ;

   2. met la somme à la puissance puissance.

   Dans un second temps, donnez à puissance la valeur 2 par défaut et faites en sorte que la fonction prenne en charge les vecteurs x présentant des valeurs manquantes NA.

   Afficher/masquer la solution

   ---

   # Associée à `<-`, la fonction function()
   # permet de définir une nouvelle fonction
   monCalcul <- function(x, puissance){
     resultat <- sum(x)^puissance
     return(resultat)
   }
   monCalcul(1:3, puissance = 2)
   ## [1] 36
   
   # La syntaxe suivante est équivalente : 
   monCalcul <- function(x, puissance) sum(x)^puissance
   monCalcul(1:3, puissance = 2)
   ## [1] 36
   
   # Pour ajouter un argument par défaut, il suffit 
   # de le spécifier dans la parenthèse de function()
   monCalcul <- function(x, puissance = 2) sum(x)^puissance
   monCalcul(1:3)
   ## [1] 36
   
   # L'argument na.rm = TRUE de la fonction sum()
   # permet d'exclure automatiquement les valeurs manquantes.
   monCalcul <- function(x, puissance = 2) sum(x, na.rm = TRUE)^puissance
   monCalcul(c(1:3, NA))
   ## [1] 36

   ---

    

3. Modifier la fonction pour ajouter une étape de vérification du type de x : prévoyez un message d’erreur si x est de type character ou factor.

   Afficher/masquer la solution

   ---

   # Les fonctions is.character() et is.factor() permettent
   # de tester le type de x. 
   monCalcul <- function(x, puissance = 2){
     if(is.character(x) || is.factor(x)) stop("x est de type caractère ou factor.")
     sum(x, na.rm = TRUE)^puissance
   }
   monCalcul(c("a", "c"))
   ## Error in monCalcul(c("a", "c")): x est de type caractère ou factor.

   ---

    

4. Quelle est la valeur de l’objet T ? Comment expliquez-vous que cet objet soit défini alors que vous ne l’avez pas vous-même créé (vous pouvez utiliser la fonction getAnywhere() pour répondre) ?

   Afficher/masquer la solution

   ---

   # L'objet T a par défaut la valeur TRUE
   T
   ## [1] TRUE
   
   # T est défini dans le package base de R comme un alias de TRUE
   getAnywhere(T)
   ## A single object matching 'T' was found
   ## It was found in the following places
   ##   package:base
   ##   namespace:base
   ## with value
   ## 
   ## [1] TRUE
   base::T
   ## [1] TRUE

   ---

    

5. Soumettez le code T <- 3. Que vaut désormais l’objet T ? Pourquoi R n’accède-t-il plus à la valeur stockée par défaut (vous pouvez utiliser la fonction search() pour répondre) ? Comment accéder désormais à la valeur par défaut ? Que retenez-vous quant à l’utilisation de T et de F en lieu et place de TRUE et FALSE dans un code ?

   Afficher/masquer la solution

   ---

   T <- 3
   T 
   ## [1] 3
   # R accède en premier lieu à l'environnement global, puis
   # aux packages dans l'ordre de search()
   search()
   ##  [1] ".GlobalEnv"             "package:foreign"       
   ##  [3] "package:bindrcpp"       "package:scales"        
   ##  [5] "package:Rcpp"           "package:Matrix"        
   ##  [7] "package:sqldf"          "package:RSQLite"       
   ##  [9] "package:gsubfn"         "package:proto"         
   ## [11] "package:microbenchmark" "package:ggplot2"       
   ## [13] "package:pryr"           "package:data.table"    
   ## [15] "package:nycflights13"   "package:dplyr"         
   ## [17] "package:MASS"           "package:ctv"           
   ## [19] "package:formatR"        "package:stats"         
   ## [21] "package:graphics"       "package:grDevices"     
   ## [23] "package:utils"          "package:datasets"      
   ## [25] "package:methods"        "Autoloads"             
   ## [27] "package:base"
   # Le package base est situé en dernière position.
   
   # Pour accéder explicitement à l'élément du package base,
   # on utilise ::
   base::T
   ## [1] TRUE
   
   # Le fait que T et F ne soit que des alias pour TRUE et FALSE
   # et donc que leur valeur soit modifiable invite à la plus
   # grande prudence dans leur utilisation. Elle est en fait
   # à proscrire dans la réalisation de projets importants 
   # (par exemple l'écriture d'un package).

   ---

    

 

Cas pratique 2 apply() : Appliquer une fonction selon les dimensions d’une matrice

1. Créez une matrice e1 de 3 lignes et de 5 colonnes en utilisant la fonction runif(). Ses valeurs sont-elles identiques si vous la générez une seconde fois ? Comment faire pour que cela soit le cas ?

   Afficher/masquer la solution

   ---

   # Par défaut,  l'utilisation sucessive de la fonction
   # runif() ne conduit pas aux mêmes résultats
   e1 <- matrix(runif(15), ncol = 5)
   e1
   ##           [,1]       [,2]      [,3]      [,4]      [,5]
   ## [1,] 0.9517595 0.09967682 0.1995403 0.6697140 0.6814065
   ## [2,] 0.5471749 0.97088164 0.2908252 0.4237895 0.6609128
   ## [3,] 0.9269703 0.60405585 0.7520751 0.9651914 0.7546232
   e1b <- matrix(runif(15), ncol = 5)
   e1b
   ##           [,1]         [,2]      [,3]       [,4]       [,5]
   ## [1,] 0.8563281 0.0129960224 0.6935136 0.87513154 0.24786063
   ## [2,] 0.7438388 0.2310905347 0.8141819 0.55749948 0.06994223
   ## [3,] 0.5767685 0.0006464222 0.1505970 0.07619284 0.53948312
   identical(e1, e1b)
   ## [1] FALSE
   
   # La fonction set.seed() permet en revanche d'initialiser
   # le générateur de nombres pseudo-aléatoires de R. 
   # Utilisé après la fonction set.seed(), deux fonctions
   # runif() donneront toujours le même résultat.
   set.seed(2016)
   e1 <- matrix(runif(15), ncol = 5)
   e1
   ##           [,1]      [,2]        [,3]       [,4]      [,5]
   ## [1,] 0.1801636 0.1335746 0.616631283 0.05345881 0.2225800
   ## [2,] 0.1429437 0.4775025 0.890547575 0.38868792 0.8765744
   ## [3,] 0.8416465 0.1212584 0.002623455 0.27295359 0.2466688
   set.seed(2016)
   e1b <- matrix(runif(15), ncol = 5)
   e1b
   ##           [,1]      [,2]        [,3]       [,4]      [,5]
   ## [1,] 0.1801636 0.1335746 0.616631283 0.05345881 0.2225800
   ## [2,] 0.1429437 0.4775025 0.890547575 0.38868792 0.8765744
   ## [3,] 0.8416465 0.1212584 0.002623455 0.27295359 0.2466688
   identical(e1, e1b)
   ## [1] TRUE

   ---

    

2. Utilisez la fonction apply() pour calculer la somme des termes de chaque ligne de e1. Comment auriez-vous pu faire autrement ? Utilisez la fonction microbenchmark() pour comparer ces deux solutions.

   Afficher/masquer la solution

   ---

   # La fonction apply() permet d'appliquer
   # une même fonction aux lignes ou aux colonnes
   # d'une matrice. 
   apply(e1, 1, sum)
   ## [1] 1.206408 2.776256 1.485151
   # Le second argument indique la dimension selon
   # laquelle appliquer la fonction (1 pour les lignes, 
   # 2 pour les colonnes). 
   
   # Cette opération est en fait nativement implémentée
   # via la fonction rowSums()
   rowSums(e1)
   ## [1] 1.206408 2.776256 1.485151
   
   # La fonction microbenchmark() permet de comparer
   # systématiquement ces deux stratégies (apply() ou
   # rowSums())
   microbenchmark(times = 1e4
     , apply = apply(e1, 1, sum)
     , rowSums = rowSums(e1)
   )
   ## Unit: microseconds
   ##     expr    min     lq    mean  median      uq      max neval
   ##    apply 21.308 24.342 28.2082 25.3315 26.6960 2252.306 10000
   ##  rowSums 16.411 19.057 21.6968 20.0705 21.3445  169.443 10000

   ---

    

3. Utilisez la fonction apply() pour centrer-réduire toutes les colonnes de e1 (i.e. leur soustraire leur moyenne puis les diviser par leur écart-type).

   Afficher/masquer la solution

   ---

   # Etape 1 : fonction pour centrer-réduire une variable
   x <- runif(10)
   mean(x) # Moyenne de x
   ## [1] 0.3157225
   sd(x) # Ecart-type de x
   ## [1] 0.2328214
   centrer_reduire <- function(x) ( x - mean(x) ) / sd(x)
   z <- centrer_reduire(x)
   mean(z)
   ## [1] 5.551115e-17
   sd(z)
   ## [1] 1
   
   # Etape 2 : utilisation dans un apply()
   e2 <- apply(e1, 2, centrer_reduire)
   colMeans(e2)
   ## [1]  0.000000e+00  3.700743e-17  1.110223e-16 -1.850372e-17  7.401487e-17
   apply(e2, 2, sd)
   ## [1] 1 1 1 1 1
   
   # Tout en une seule étape :
   e3 <- apply(e1, 2, function(x) ( x - mean(x) ) / sd(x))
   identical(e2, e3)
   ## [1] TRUE

   ---

    

 

Cas pratique 3 lapply() et sapply() : Appliquer une fonction aux éléments d’un vecteur, d’une liste ou aux colonnes d’un data.frame

1. On définit le vecteur f1 par f1 <- 5:15. Utilisez la fonction sapply() pour calculer la somme cumulée des éléments de f1. Quelle(s) alternative(s) envisageriez-vous ? Utilisez la fonction microbenchmark() pour comparer ces solutions.

   Afficher/masquer la solution

   ---

   # INDICATION : commencer par définir la fonction 
   # sumfirst(x, i) qui calcule la somme des i premiers 
   # éléments de x, puis utilisez-la dans un sapply(). 
   
   f1 <- 5:15
   
   # Méthode 1 : sapply()
   # L'objectif est d'obtenir un vecteur de longueur length(f1)
   # qui pour chaque élément i renvoie la somme des éléments 1 à i
   # de f1
   
   # Pour calculer la somme des 4 premiers termes, il suffit d'appliquer
   # la fonction sum() au vecteur f1 restreint à ses 4 premiers termes
   sum(f1[1:4])
   ## [1] 26
   
   # On généralise cette idée dans un sapply() en définissant
   # la fonction à la volée (x remplace ici le 4 de l'exemple)
   sapply(1:length(f1), function(i) sum(x[1:i]))
   ##  [1] 0.5808441 1.2996963 1.7508546 1.9072832 1.9705809 1.9829998 2.4652852
   ##  [8] 2.7432785 2.9876747 3.1572245        NA
   # Note : dans cette méthode pour chaque élément du vecteur
   # on recalcule la somme, sans réutiliser celle calculée
   # pour les éléments précédents.
   
   # Remarque : on peut aussi utiliser la fonction
   # seq_along() pour créer le vecteur d'indices
   sapply(seq_along(f1), function(i) sum(x[1:i]))
   ##  [1] 0.5808441 1.2996963 1.7508546 1.9072832 1.9705809 1.9829998 2.4652852
   ##  [8] 2.7432785 2.9876747 3.1572245        NA
   
   
   
   # Méthode 2 : cumsum()
   # La fonction cumsum() effectue nativement cette opération
   cumsum(f1)
   ##  [1]   5  11  18  26  35  45  56  68  81  95 110
   
   # Comparaison avec microbenchmark()
   microbenchmark(times = 1e4
     , sapply = sapply(seq_along(f1), function(i) sum(f1[1:i]))
     , cumsum = cumsum(f1)
   )
   ## Unit: nanoseconds
   ##    expr   min    lq       mean  median    uq     max neval
   ##  sapply 27908 29722 34932.9104 31192.5 33480 5585155 10000
   ##  cumsum   260   367   549.6245   545.0   629   15181 10000

   ---

    

2. On définit la liste f2 par

   f2 <- list(
     sample.int(26, 10, replace = TRUE)
     , sample.int(26, 100, replace = TRUE)
     , sample.int(26, 1000, replace = TRUE)
   )

Utilisez les fonctions lapply() ou sapply() pour :

• retrouver la longueur de chacun des éléments de f2 ;
• extraire les 5 premiers éléments de chacun des éléments de f2 ;

• remplacer chaque élément de f2 par le vecteur de lettres (en minuscules) dont il représente les positions dans l’alphabet.

  Afficher/masquer la solution

  ---

  # Les fonctions sapply() et lapply() permettent 
  # d'appliquer systématiquement une fonction aux éléments
  # d'une liste. 
  
  # Pour connaître la longueur de chaque élément de f2, 
  # il suffit ainsi d'appliquer à chacun la fonction length()
  lapply(f2, length)
  ## [[1]]
  ## [1] 10
  ## 
  ## [[2]]
  ## [1] 100
  ## 
  ## [[3]]
  ## [1] 1000
  
  # sapply() effectue exactement le même traitement que
  # lapply() mais essaie en plus de simplifier le résultat
  # si c'est possible
  sapply(f2, length)
  ## [1]   10  100 1000
  # Ici c'est le cas : on passe d'une liste avec lapply()
  # à un vecteur avec sapply()
  
  # Pour extraire les éléments 1 à 10 de chacun des éléments
  # de f2, on utilise l'opérateur [
  lapply(f2, function(x) x[1:10])
  ## [[1]]
  ##  [1]  2 17 13 12 25 25 17 16  4 20
  ## 
  ## [[2]]
  ##  [1]  4  3  4  9 11  8  7  9 10 17
  ## 
  ## [[3]]
  ##  [1] 18  5 15  8 14  6  7 13 26  2
  # On aurait aussi pu soumettre de façon équivalente 
  lapply(f2, `[`, 1:10)
  ## [[1]]
  ##  [1]  2 17 13 12 25 25 17 16  4 20
  ## 
  ## [[2]]
  ##  [1]  4  3  4  9 11  8  7  9 10 17
  ## 
  ## [[3]]
  ##  [1] 18  5 15  8 14  6  7 13 26  2
  # car x[1:10] est équivalent à `[`(x, 1:10)
  
  # Remplacer les éléments d'un des éléments de f2
  # par les lettres dont il représente la position
  # est relativement simple : 
  letters[f2[[1]]]
  ##  [1] "b" "q" "m" "l" "y" "y" "q" "p" "d" "t"
  # On reprend donc le principe de la sous-question
  # précédente dans un lapply()
  # lapply(f2, function(x) letters[x])
  # (commenté ici pour ne pas saturer la sortie).

  ---

   

3. On définit le data.frame f3 par

   set.seed(1)
   f3 <- data.frame(
     id = letters[1:20]
     , by = rep(letters[1:5], times = 4)
     , matrix(runif(100), ncol = 5)
     , stringsAsFactors = FALSE
   )

   Utilisez les fonctions lapply() ou sapply() pour :

• déterminer le type de chacune des variables de f3;
• calculer la moyenne de toutes les variables numériques de f3;
• convertir toutes les variables de type character en facteurs.
  Afficher/masquer la solution

  ---

  # Un data.frame étant un cas particulier de liste, 
  # on peut lui appliquer les fonctions lapply() et
  # sapply() exactement comme on le ferait pour une liste.
  
  # Pour connaître le type de chacun de ses éléments, 
  # on leur applique donc systématiquement la fonction
  # typeof()
  sapply(f3, typeof)
  ##          id          by          X1          X2          X3          X4 
  ## "character" "character"    "double"    "double"    "double"    "double" 
  ##          X5 
  ##    "double"
  
  # Pour appliquer une même fonction à toutes ses variables 
  # de type numériques, on commence par les identifier
  # avec un is.numeric() dans un sapply()
  num <- sapply(f3, is.numeric)
  num
  ##    id    by    X1    X2    X3    X4    X5 
  ## FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
  # Puis on leur applique la fonction désirée
  sapply(f3[num], mean)
  ##        X1        X2        X3        X4        X5 
  ## 0.5551671 0.4738822 0.5072134 0.5665520 0.4864206
  
  # De même pour les variables de type caractère
  char <- sapply(f3, is.character)
  char
  ##    id    by    X1    X2    X3    X4    X5 
  ##  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE
  f3[char] <- lapply(f3[char], as.factor)
  str(f3)
  ## 'data.frame':    20 obs. of  7 variables:
  ##  $ id: Factor w/ 20 levels "a","b","c","d",..: 1 2 3 4 5 6 7 8 9 10 ...
  ##  $ by: Factor w/ 5 levels "a","b","c","d",..: 1 2 3 4 5 1 2 3 4 5 ...
  ##  $ X1: num  0.266 0.372 0.573 0.908 0.202 ...
  ##  $ X2: num  0.935 0.212 0.652 0.126 0.267 ...
  ##  $ X3: num  0.821 0.647 0.783 0.553 0.53 ...
  ##  $ X4: num  0.913 0.294 0.459 0.332 0.651 ...
  ##  $ X5: num  0.435 0.713 0.4 0.325 0.757 ...

  ---

   

 

Cas pratique 4 tapply() : Appliquer une fonction selon les modalités d’un vecteur

1. On définit le vecteur g1 <- sample(20) et le vecteur g2 <- rep(c("H", "F"), times = 10). Utilisez tapply() pour calculer la moyenne de g1 selon les groupes définis par g2.
   Afficher/masquer la solution

   ---

   set.seed(2)
   g1 <- sample(20)
   g2 <- rep(c("H", "F"), times = 10)
   
   # La fonction tapply() permet d'appliquer une fonction
   # à un vecteur selon les modalités d'un autre vecteur
   tapply(g1, g2, mean)
   ##    F    H 
   ## 10.7 10.3

   ---

    

2. On repart du data.frame f3 du cas pratique précédent. Utilisez tapply() pour calculer le total de la variable X1 selon les modalités de la variable by.
   Afficher/masquer la solution

   ---

   # Il suffit d'utiliser tapply() sur les variables de f3
   tapply(f3$X1, f3$by, sum)
   ##        a        b        c        d        e 
   ## 1.867572 2.210974 2.912580 2.301461 1.810755

   ---

    

3. Combinez la fonction split() avec sapply() pour obtenir le même résultat. Comment calculeriez-vous le total de toutes les variables numériques de f3 selon les modalités de la variable by ?

   Afficher/masquer la solution

   ---

   # La fonction split() éclate en une liste un vecteur 
   # ou un data.frame selon les modalités d'une ou plusieurs
   # variables.
   f3[,c("by", "X1")]
   ##    by         X1
   ## 1   a 0.26550866
   ## 2   b 0.37212390
   ## 3   c 0.57285336
   ## 4   d 0.90820779
   ## 5   e 0.20168193
   ## 6   a 0.89838968
   ## 7   b 0.94467527
   ## 8   c 0.66079779
   ## 9   d 0.62911404
   ## 10  e 0.06178627
   ## 11  a 0.20597457
   ## 12  b 0.17655675
   ## 13  c 0.68702285
   ## 14  d 0.38410372
   ## 15  e 0.76984142
   ## 16  a 0.49769924
   ## 17  b 0.71761851
   ## 18  c 0.99190609
   ## 19  d 0.38003518
   ## 20  e 0.77744522
   split(f3$X1, f3$by)
   ## $a
   ## [1] 0.2655087 0.8983897 0.2059746 0.4976992
   ## 
   ## $b
   ## [1] 0.3721239 0.9446753 0.1765568 0.7176185
   ## 
   ## $c
   ## [1] 0.5728534 0.6607978 0.6870228 0.9919061
   ## 
   ## $d
   ## [1] 0.9082078 0.6291140 0.3841037 0.3800352
   ## 
   ## $e
   ## [1] 0.20168193 0.06178627 0.76984142 0.77744522
   
   # Ce faisant, on peut avec une fonction lapply()
   # ou sapply() reproduire le comportent de la fonction
   # tapply()
   sapply(split(f3$X1, f3$by), sum)
   ##        a        b        c        d        e 
   ## 1.867572 2.210974 2.912580 2.301461 1.810755
   
   # L'avantage de cette technique est qu'elle
   # permet d'appliquer le même type de traitement
   # non à une seule variable mais à plusieurs
   # d'un seul coup
   sapply(split(f3[num], f3$by), function(x){
     sapply(x, sum)
   })
   ##           a        b        c        d        e
   ## X1 1.867572 2.210974 2.912580 2.301461 1.810755
   ## X2 2.471366 1.619339 1.635547 1.905174 1.846217
   ## X3 2.187388 1.847873 2.216894 2.192152 1.699960
   ## X4 2.402164 2.475928 1.962049 1.527737 2.963161
   ## X5 1.674290 1.937845 1.574059 2.257980 2.284239
   # ... ou de façon équivalente (mais un peu difficile
   # à relire !)
   sapply(split(f3[num], f3$by), sapply, sum)
   ##           a        b        c        d        e
   ## X1 1.867572 2.210974 2.912580 2.301461 1.810755
   ## X2 2.471366 1.619339 1.635547 1.905174 1.846217
   ## X3 2.187388 1.847873 2.216894 2.192152 1.699960
   ## X4 2.402164 2.475928 1.962049 1.527737 2.963161
   ## X5 1.674290 1.937845 1.574059 2.257980 2.284239
   # ... ou plus efficacement avec colSums()
   sapply(split(f3[num], f3$by), colSums)
   ##           a        b        c        d        e
   ## X1 1.867572 2.210974 2.912580 2.301461 1.810755
   ## X2 2.471366 1.619339 1.635547 1.905174 1.846217
   ## X3 2.187388 1.847873 2.216894 2.192152 1.699960
   ## X4 2.402164 2.475928 1.962049 1.527737 2.963161
   ## X5 1.674290 1.937845 1.574059 2.257980 2.284239

   ---

    

 

Cas pratique 5 do.call() : Appliquer une fonction simultanément à l’ensemble des éléments d’une liste

De nombreuses fonctions peuvent porter sur un nombre indéterminé d’éléments : c(), sum(), rbind(), etc. Pour les appliquer à l’ensemble des éléments d’une liste sans avoir à tous les écrire un à un, il suffit d’utiliser do.call().

1. On définit la liste h1 <- list(1:5, 6:10, 11:15). Que se passe-t-il si vous soumettez sum(h1) ? Utilisez do.call() pour sommer l’ensemble des éléments de h1. Comparez avec le résultat de sapply(h1, sum).

   Afficher/masquer la solution

   ---

   h1 <- list(1:5, 6:10, 11:15)
   
   # La fonction sum() ne peut pas porter sur une liste
   sum(h1)
   ## Error in sum(h1): 'type' (list) de l'argument incorrect
   
   # En revanche, il est possible d'appliquer la fonction
   # sum() à l'ensemble des éléments de h1
   sum(h1[[1]], h1[[2]], h1[[3]])
   ## [1] 120
   
   # Pour automatiser cette expression, on peut utiliser
   # do.call()
   do.call(sum, h1)
   ## [1] 120
   
   # A l'inverse, sapply(h1, sum) applique la fonction
   # sum() à chaque élément de h1 : 
   sapply(h1, sum)
   ## [1] 15 40 65

   ---

    

2. Réunissez les éléments de h1 en un seul vecteur avec la fonction base::c(). Comparez avec le résultat de lapply(h1, base::c)

   Afficher/masquer la solution

   ---

   # La fonction base::c() permet de concaténer
   # les vecteurs qui constituent h1 
   c(h1[[1]], h1[[2]], h1[[3]])
   ##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
   
   # Pour automatiser cette expression, on peut utiliser
   # do.call()
   do.call(base::c, h1)
   ##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
   # Note : on utilise ici base::c pour bien indiquer que
   # c'est la fonction c() du package base que l'on 
   # souhaite utiliser. En effet, si l'on avait défini
   # on objet c dans l'environnement global il y aurait
   # eu un conflit de noms.
   
   # lapply(h1, base::c) applique la fonction base::c()
   # à chaque élément de h1, ce qui ne change rien
   lapply(h1, base::c)
   ## [[1]]
   ## [1] 1 2 3 4 5
   ## 
   ## [[2]]
   ## [1]  6  7  8  9 10
   ## 
   ## [[3]]
   ## [1] 11 12 13 14 15

   ---

    

3. On définit la liste de matrices h2

   h2 <- list(
     matrix(1:6, nrow = 2)
     , matrix(7:12, nrow = 2)
     , matrix(13:18, nrow = 2)
   )

   Utilisez do.call() avec les fonctions rbind() et cbind() pour concaténer l’ensemble des éléments de h2 en ligne ou en colonne respectivement.
   Afficher/masquer la solution

   ---

   # Pour concaténer les éléments de h2 en ligne, 
   # on utilise rbind()
   rbind(h2[[1]], h2[[2]], h2[[3]])
   ##      [,1] [,2] [,3]
   ## [1,]    1    3    5
   ## [2,]    2    4    6
   ## [3,]    7    9   11
   ## [4,]    8   10   12
   ## [5,]   13   15   17
   ## [6,]   14   16   18
   
   # Pour automatiser cette expression, on peut utiliser
   # do.call()
   do.call(rbind, h2)
   ##      [,1] [,2] [,3]
   ## [1,]    1    3    5
   ## [2,]    2    4    6
   ## [3,]    7    9   11
   ## [4,]    8   10   12
   ## [5,]   13   15   17
   ## [6,]   14   16   18
   
   # De même en colonne avec cbind()
   do.call(cbind, h2)
   ##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
   ## [1,]    1    3    5    7    9   11   13   15   17
   ## [2,]    2    4    6    8   10   12   14   16   18

   ---

    

 

Cas pratique 6 Reduce() : Appliquer une fonction successivement à l’ensemble des éléments d’une liste

De nombreuses fonctions portent sur deux éléments précisément : opérations arithmétiques, merge(), etc. Pour les appliquer successivement à l’ensemble des éléments d’une liste, il suffit d’utiliser Reduce().

1. On repart de la liste h1 du cas pratique précédent. Observez le résultat de Reduce(`+`, h1) : comment le comprenez-vous ? Comparez avec sapply(h1, sum).

   Afficher/masquer la solution

   ---

   # Reduce(f, x) applique la fonction f() aux éléments de x
   # 2 à 2 et successivement. 
   Reduce(`+`, h1)
   ## [1] 18 21 24 27 30
   # Ici, cette expression est équivalente à 
   ( h1[[1]] + h1[[2]] ) + h1[[3]]
   ## [1] 18 21 24 27 30
   # Le vecteur obtenu est la somme terme à terme de 
   # tous les éléments qui composent h1
   
   # A l'inverse, sapply(h1, sum) applique la fonction
   # sum() à chaque élément de h1 : 
   sapply(h1, sum)
   ## [1] 15 40 65

   ---

    

2. On définit la liste de data.frame i1

   i1 <- list(
     data.frame(id = letters[1:4], var1 = 1:4, stringsAsFactors = FALSE)
     , data.frame(id = letters[2:5], var2 = 5:8, stringsAsFactors = FALSE)
     , data.frame(id = letters[3:6], var3 = 9:12, stringsAsFactors = FALSE)
   )

   Utilisez Reduce() pour fusionner l’ensemble des éléments de i1 selon la variable id. Comment ajusteriez-vous ce code pour pouvoir utiliser l’option all = TRUE de la fonction merge() ?
   Afficher/masquer la solution

   ---

   # La fonction Reduce() est particulièrement utile 
   # pour fusionner de nombreux data.frame selon la
   # même variable. Par défaut, la fonction merge()
   # fusionne les data.frame sur les variables qu'ils
   # ont en commun.
   Reduce(merge, i1)
   ##   id var1 var2 var3
   ## 1  c    3    6    9
   ## 2  d    4    7   10
   # Ce code fournit le résultat désiré. 
   
   # Pour spécifier explicitement la variable de fusion, 
   # il faut redéfinir à la volée la fonction à appliquer 
   # dans le Reduce() (comme dans un *apply) : 
   Reduce(function(x, y) merge(x, y, by = "id"), i1)
   ##   id var1 var2 var3
   ## 1  c    3    6    9
   ## 2  d    4    7   10
   
   # De même pour ajouter d'autres options à la fonction 
   # merge()
   Reduce(function(x, y) merge(x, y, by = "id", all = TRUE), i1)
   ##   id var1 var2 var3
   ## 1  a    1   NA   NA
   ## 2  b    2    5   NA
   ## 3  c    3    6    9
   ## 4  d    4    7   10
   ## 5  e   NA    8   11
   ## 6  f   NA   NA   12

   ---

    

Travailler efficacement sur des données avec base R

L’ensemble des cas pratiques qui suivent portent sur les données de l’enquête Emploi en continu stockées dans le fichier eect4.rds. La fonction readRDS() permet de le charger en mémoire :

setwd("Y:/Documentation/R/R_perfectionnement/donnees")
eec <- readRDS("eect4.rds")
str(eec)

## 'data.frame':    34913 obs. of  19 variables:
##  $ IDENT    : chr  "G0A56JP6" "G0A56JP6" "G0A56JR6" "G0A56JS6" ...
##  $ TRIM     : chr  "4" "4" "4" "4" ...
##  $ NOI      : chr  "01" "02" "01" "01" ...
##  $ REG      : chr  "11" "11" "11" "11" ...
##  $ AGE      : chr  "66" "29" "27" "29" ...
##  $ SEXE     : chr  "2" "1" "2" "2" ...
##  $ CSE      : chr  "56" "81" "38" "37" ...
##  $ DIP11    : chr  "71" "42" "10" "11" ...
##  $ ACTEU    : chr  "1" "2" "1" "1" ...
##  $ SALRED   : int  596 NA 2700 2666 11967 NA 2000 2800 2333 3500 ...
##  $ STC      : chr  "2" NA "2" "2" ...
##  $ TAM1D    : chr  NA NA NA NA ...
##  $ AIDREF   : chr  NA "5" NA NA ...
##  $ TPP      : chr  "1" NA "1" "1" ...
##  $ NBAGENF  : chr  "0" "0" "0" "0" ...
##  $ DUHAB    : chr  "7" NA "7" "7" ...
##  $ PUB3FP   : chr  "4" NA "4" "4" ...
##  $ NAIA     : chr  "1946" "1983" "1985" "1983" ...
##  $ EXTRI1613: num  1777 1777 2045 1898 1754 ...

 

Cas pratique 7 Sélection d’observations

1. Utilisez l’opérateur [ pour sélectionner l’individu appartenant au ménage dont l’IDENT est "GFO5NVUE" et dont le numéro d’ordre NOI est "02".

   Afficher/masquer la solution

   ---

   # Il suffit de créer le vecteur logique correspondant
   # et de l'utiliser das `[`
   eec[eec$IDENT == "GFO5NVUE" & eec$NOI == "02", ]  
   ##           IDENT TRIM NOI REG AGE SEXE CSE DIP11 ACTEU SALRED STC TAM1D
   ## 360720 GFO5NVUE    4  02  42  56    1  63    50     1   2004   2  <NA>
   ##        AIDREF TPP NBAGENF DUHAB PUB3FP NAIA EXTRI1613
   ## 360720   <NA>   1       0     6      4 1956  990.2715

   ---

    

2. Cherchez de la documentation sur la fonction subset(). Comparez ses performances avec celles de l’opérateur [ à l’aide de la fonction microbenchmark().

   Afficher/masquer la solution

   ---

   # subset() permet d'évaluer une clause logique
   # sans avoir à répéter le nom du data.frame
   # d'origine
   subset(eec, IDENT == "GFO5NVUE" & NOI == "02")
   ##           IDENT TRIM NOI REG AGE SEXE CSE DIP11 ACTEU SALRED STC TAM1D
   ## 360720 GFO5NVUE    4  02  42  56    1  63    50     1   2004   2  <NA>
   ##        AIDREF TPP NBAGENF DUHAB PUB3FP NAIA EXTRI1613
   ## 360720   <NA>   1       0     6      4 1956  990.2715
   
   # Comparaison des performances
   microbenchmark(times = 1e2
     , "[" = eec[eec$IDENT == "GFO5NVUE" & eec$NOI == "02", ]  
     , subset = subset(eec, IDENT == "GFO5NVUE" & NOI == "02")
   )
   ## Unit: milliseconds
   ##    expr      min       lq     mean   median       uq      max neval
   ##       [ 1.323394 1.386014 1.640976 1.426236 1.501453 7.292565   100
   ##  subset 1.497029 1.555517 1.919099 1.613505 1.712313 7.904183   100

   ---

    

3. Concaténez les valeurs des variables IDENT et NOI et utilisez le résultat comme noms de ligne. Retrouvez alors l’individu de la question a à l’aide de son nom de ligne. Comparez les performances de cette méthode avec celles de l’opérateur [.

   Afficher/masquer la solution

   ---

   # La foncton paste0() permet de concaténer des chaînes
   # de caractères sans utiliser de délimiteur
   row.names(eec) <- paste0(eec$IDENT, eec$NOI)
   
   # Sélection de l'individu par son nom de ligne
   eec["GFO5NVUE02", ]
   ##               IDENT TRIM NOI REG AGE SEXE CSE DIP11 ACTEU SALRED STC TAM1D
   ## GFO5NVUE02 GFO5NVUE    4  02  42  56    1  63    50     1   2004   2  <NA>
   ##            AIDREF TPP NBAGENF DUHAB PUB3FP NAIA EXTRI1613
   ## GFO5NVUE02   <NA>   1       0     6      4 1956  990.2715
   
   # Comparaison des performances
   microbenchmark(times = 1e2
     , "[" = eec[eec$IDENT == "GFO5NVUE" & eec$NOI == "02", ]  
     , names = eec["GFO5NVUE02", ]
   )
   ## Unit: microseconds
   ##   expr      min       lq     mean   median       uq      max neval
   ##      [ 1320.147 1428.255 1642.419 1500.146 1553.541 9371.659   100
   ##  names  896.914 1028.111 1217.348 1146.683 1232.324 8138.571   100

   ---

    

 

Cas pratique 8 Création de variables

On souhaite recoder la variable AGE en trois modalités : 15-30 ans, 31-60 ans et plus de 60 ans. On part du code suivant :

for(i in 1:nrow(eec)) eec$trage1[i] <- if(as.numeric(eec$AGE[i]) < 31) "15-30" else if(as.numeric(eec$AGE[i]) < 61) "31-60" else "61 et +"

1. Mesurez le temps d’exécution du code proposé. Que pensez-vous de cette syntaxe ?

   Afficher/masquer la solution

   ---

   # Mesure du temps d'exécution de cette première version
   microbenchmark(times = 1
     , v1 = {
   for(i in 1:nrow(eec)) eec$trage1[i] <- if(as.numeric(eec$AGE[i]) < 31) "15-30" else if(as.numeric(eec$AGE[i]) < 61) "31-60" else "61 et +"
     }
   )
   ## Unit: seconds
   ##  expr      min       lq     mean   median       uq      max neval
   ##    v1 7.826031 7.826031 7.826031 7.826031 7.826031 7.826031     1
   # Remarques : 
   # - on ne fait qu'une seule itération car le temps d'exécution 
   # est très long (plusieurs secondes !).
   # - le as.numeric() n'est pas indispensable mais accélère sensiblement
   # le traitement en faisant en sorte que la comparaison par < s'effectue 
   # entre vecteurs de type numérique et pas entre vecteurs de type caractère.
   
   # Cette syntaxe présente un énorme problème : elle utilise
   # une boucle là où des opérations vectorisées beaucoup 
   # plus rapides sont disponibles.

   ---

    

2. On propose une seconde version du code :

   eec$trage2 <- ifelse(as.numeric(eec$AGE) < 31, "15-30", ifelse(
     as.numeric(eec$AGE) < 61, "31-60", "61 et +"
   ))

   Vérifiez que le résultat est identique à la proposition initiale (par exemple avec identical()ou all.equal()) et mesurez le temps temps d’exécution de cette deuxième version. Êtes-vous surpris du gain de performances ?

   Afficher/masquer la solution

   ---

   eec$trage2 <- ifelse(as.numeric(eec$AGE) < 31, "15-30", ifelse(
     as.numeric(eec$AGE) < 61, "31-60", "61 et +"
   ))
   
   # Vérification de la correspondance 
   # avec la première version
   all.equal(eec$trage1, eec$trage2)
   ## [1] TRUE
   
   # Comparaison des performances
   microbenchmark(times = 1e1
     , v2 = {
   eec$trage2 <- ifelse(as.numeric(eec$AGE) < 31, "15-30", ifelse(
     as.numeric(eec$AGE) < 61, "31-60", "61 et +"
   ))
     }
   )
   ## Unit: milliseconds
   ##  expr      min       lq     mean   median       uq      max neval
   ##    v2 22.00497 22.06697 22.09491 22.10118 22.12497 22.20733    10
   
   # On n'est pas surpris du gain de performance : en R, 
   # les boucles sont beaucoup moins efficaces que les 
   # opérations vectorisées (qui reposent sur des 
   # boucles dans des langages de plus bas niveau).

   ---

    

3. Comment recoderiez-vous la proposition précédente pour ne plus faire appel à la fonction ifelse() ? Cela est-il susceptible d’améliorer les performances ? Mettez en oeuvre cette solution et mesurez-en les performances.

   Afficher/masquer la solution

   ---

   # Il est possible de se passer de ifelse() en utilisant
   # l'opérateur [ pour effectuer des remplacements successifs. 
   eec$trage3 <- "15-30"
   eec$trage3[as.numeric(eec$AGE) > 30 & as.numeric(eec$AGE) < 61] <- "31-60"
   eec$trage3[as.numeric(eec$AGE) > 60] <- "61 et +"
   
   # Vérification de la correspondance 
   # avec la première version
   all.equal(eec$trage1, eec$trage3)
   ## [1] TRUE
   
   # Comparaison des performances
   microbenchmark(times = 1e1
     , v2 = {
   eec$trage2 <- ifelse(as.numeric(eec$AGE) < 31, "15-30", ifelse(
     as.numeric(eec$AGE) < 61, "31-60", "61 et +"
   ))
     }
     , v3 = {
   eec$trage3 <- "15-30"
   eec$trage3[as.numeric(eec$AGE) > 30 & as.numeric(eec$AGE) < 61] <- "31-60"
   eec$trage3[as.numeric(eec$AGE) > 60] <- "61 et +"
     }
   )
   ## Unit: milliseconds
   ##  expr       min        lq     mean    median        uq      max neval
   ##    v2 21.932762 21.981801 22.77927 22.029922 22.285772 29.08888    10
   ##    v3  9.956607  9.970045 10.00405  9.982753  9.997636 10.18163    10

   ---

    

4. Combien de fois appelez-vous la fonction as.numeric() dans le code qui précède ? Proposez une dernière version qui minimise les opérations effectuées par R et synthétisez le gain en termes de performances.

   Afficher/masquer la solution

   ---

   # La fonction as.numeric() est appelée trois fois dans
   # le code qui précède. Pour gagner du temps, on
   # crée la variable t temporaire que l'on utilise
   # en lieu et place de as.numeric(eec$AGE)
   t <- as.numeric(eec$AGE)
   eec$trage4 <- "15-30"
   eec$trage4[t > 30 & t < 61] <- "31-60"
   eec$trage4[t > 60] <- "61 et +"
   
   # Vérification de la correspondance 
   # avec la première version
   all.equal(eec$trage1, eec$trage4)
   ## [1] TRUE
   
   
   # Comparaison des performances
   microbenchmark(times = 10
     , v2 = {
   eec$trage2 <- ifelse(as.numeric(eec$AGE) < 31, "15-30", ifelse(
     as.numeric(eec$AGE) < 61, "31-60", "61 et +"
   ))
     }
     , v3 = {
   eec$trage3 <- "15-30"
   eec$trage3[as.numeric(eec$AGE) > 30 & as.numeric(eec$AGE) < 61] <- "31-60"
   eec$trage3[as.numeric(eec$AGE) > 60] <- "61 et +"
     }
     , v4 = {
   t <- as.numeric(eec$AGE)
   eec$trage4 <- "15-30"
   eec$trage4[t > 30 & t < 61] <- "31-60"
   eec$trage4[t > 60] <- "61 et +"
     }
   )
   ## Unit: milliseconds
   ##  expr       min        lq     mean    median        uq       max neval
   ##    v2 21.977126 21.987522 22.88556 22.032780 22.170315 30.435724    10
   ##    v3  9.960993  9.981585 10.03640 10.009300 10.071748 10.224556    10
   ##    v4  4.123373  4.126897  4.17540  4.156786  4.236318  4.268828    10
   
   # Le passage de la version 1 à la version 4 du code s'est
   # donc traduit par un gain en termes de performances
   # de l'ordre d'un facteur 1 000.

   ---

    

 

Cas pratique 9 Agrégation par groupes : salaire moyen par région

On cherche à calculer le plus efficacement possible le salaire moyen (variable SALRED) par région (variable REG), d’abord sans pondérer puis en pondérant par le poids de sondage EXTRI1613.

1. Comparez les performances des fonctions aggregate(), by(), sapply() (avec split()) et tapply() pour le calcul du salaire moyen non-pondéré par région.

   Afficher/masquer la solution

   ---

   # INDICATION : consultez l'aide de chacune de ces fonctions
   # pour vous approprier leur syntaxe.
   
   # On adapte à chaque fonction la syntaxe à mettre en oeuvre
   # et on compare leurs performances
   microbenchmark(times = 10
     , aggregate = aggregate(eec$SALRED, list(eec$REG), mean, na.rm = TRUE)
     , by = by(eec$SALRED, eec$REG, mean, na.rm = TRUE)
     , sapply = sapply(split(eec$SALRED, eec$REG), mean, na.rm = TRUE)
     , tapply = tapply(eec$SALRED, eec$REG, mean, na.rm = TRUE)
   )
   ## Unit: milliseconds
   ##       expr       min        lq      mean    median        uq       max
   ##  aggregate 16.794650 16.976295 18.175801 17.120999 17.802278 26.199475
   ##         by  3.089860  3.210694  3.315335  3.302460  3.331780  3.721400
   ##     sapply  2.020594  2.051265  2.154060  2.120222  2.165440  2.636258
   ##     tapply  2.101470  2.161448  2.180465  2.172873  2.208426  2.299784
   ##  neval
   ##     10
   ##     10
   ##     10
   ##     10
   # En règle générale sapply() et tapply() conduisent aux meilleures
   # performances et sont proches l'une de l'autre. 

   ---

    

2. Quelles pistes envisageriez-vous pour améliorer encore les performances dans ce type de situation ?
   Afficher/masquer la solution

   ---

   # Plusieurs pistes sont envisageables : 
   # - utiliser des manipulations purement vectorielles ;
   # - utiliser des manipulations matricielles sur des matrices lacunaires avec le package `Matrix`;
   # - paralléliser l'exécution avec le package `parallel`; 
   # - coder la fonction en C++.

   ---

    

3. Utilisez la fonction sapply() pour calculer le salaire moyen pondéré (par le poids de sondage EXTRI1613) par région. Y parvenez-vous également avec tapply() ?

   Afficher/masquer la solution

   ---

   # INDICATION 1 : appliquez la fonction split() à un objet
   # contenant le salaire et le poids de sondage pour 
   # pouvoir utiliser les deux variables dans le *apply(). 
   # INDICATION 2 : vous pouvez calculer la moyenne pondérée
   # "à la main" ou utiliser la fonction weighted.mean()
   
   # Par rapport à la question a., la principale différence
   # vient du fait que l'on a besoin de plusieurs éléments
   # dans chaque bloc éclaté par la fonction split() : 
   # le salaire d'une part, le poids de sondage d'autre part.
   
   # Avec la fonction weighted.mean()
   sapply(
     split(eec[c("SALRED", "EXTRI1613")], eec$REG)
     , function(x) weighted.mean(x$SALRED, x$EXTRI1613, na.rm = TRUE)
   )
   
   # Manuellement en pensant bien à exclure les poids de sondage
   # des individus pour lesquels SALRED est NA
   sapply(
     split(eec[c("SALRED", "EXTRI1613")], eec$REG)
     , function(x) sum(x$SALRED * x$EXTRI1613, na.rm = TRUE) / sum((!is.na(x$SALRED)) * x$EXTRI1613, na.rm = TRUE)
   )

   ---

    

 

Cas pratique 10 Fusion de tables : nombre d’individus au chômage par ménage

L’objectif de ce cas pratique est de créer, dans la table eec, une variable indiquant pour chaque individu le nombre d’individus au chômage dans son ménage. La position sur le marché du travail est codée par la variable ACTEU (ACTEU == "2" correspond au chômage) et l’identifiant du ménage est la variable IDENT (les individus d’un même ménage ont la même valeur pour la variable IDENT).

1. Utilisez la fonction tapply() pour déterminer le nombre de personnes au chômage dans chaque ménage et stockez cette information dans un objet appelé nbcho. Quelles sont ses caractéristiques ?

   Afficher/masquer la solution

   ---

   # Il suffit d'appliquer la fonction sum() par ménage 
   # à l'indicatrice de chômage eec$ACTEU == 2
   nbcho <- tapply(eec$ACTEU == "2", eec$IDENT, sum, na.rm = TRUE)
   
   # Caractéristiques de nbcho
   str(nbcho)
   ##  int [1:18903(1d)] 1 0 0 0 0 0 0 0 0 0 ...
   ##  - attr(*, "dimnames")=List of 1
   ##   ..$ : chr [1:18903] "G0A56JP6" "G0A56JR6" "G0A56JS6" "G0A56JT6" ...
   # nbcho est un vecteur nommé dont les noms sont 
   # les identifiants de ménage.

   ---

    

2. Utilisez la fonction merge() pour refusionner le résultat de la question précédente avec la table eec et créer la variable nbcho.

   Afficher/masquer la solution

   ---

   # Etape 1 : constituer un data.frame à partir de nbcho 
   # avec IDENT comme identifiant
   nbcho_df <- data.frame(IDENT = names(nbcho), nbcho1 = nbcho)
   
   # Etape 2 : fusionner eec et nbcho_df par IDENT
   eec <- merge(eec, nbcho_df, by = "IDENT")

   ---

    

3. Utilisez habilement les noms de vecteur et l’opérateur [ pour reproduire plus efficacement le résultat de la question b (toujours en repartant de nbcho). Vérifiez que la variable créée est bien identique.

   Afficher/masquer la solution

   ---

   # INDICATION : essayer de sélectionner par leur nom
   # les valeurs de nbcho correspondant aux 10 premières observations
   # de eec
   
   # Le point essentiel est de remarquer que l'on peut
   # utiliser les noms pour réarranger les valeurs
   # de nbcho de sorte à ce qu'elles correspondent
   # à l'ordre de eec
   
   # Identifiant des 10 premières observations de eec
   eec$IDENT[1:10]
   ##  [1] "G0A56JP6" "G0A56JP6" "G0A56JR6" "G0A56JS6" "G0A56JT6" "G0A56JU6"
   ##  [7] "G0A56JV6" "G0A56JV6" "G0A56JW6" "G0A56JX6"
   
   # Pour extraire les valeurs de nbcho correspondantes,
   # il suffit d'exploiter le fait que nbcho soit 
   # nommé à l'aide des identifiants de ménage.
   
   # Par exemple : 
   # - Valeur de nbcho pour le premier ménage de eec
   nbcho["G0A56JP6"] 
   ## G0A56JP6 
   ##        1
   # - Valeur de nbcho pour les quatre premiers ménages
   # de eec
   nbcho[c("G0A56JP6", "G0A56JP6", "G0A56JR6", "G0A56JS6")]
   ## G0A56JP6 G0A56JP6 G0A56JR6 G0A56JS6 
   ##        1        1        0        0
   # - Valeur de nbcho pour les 10 premiers ménages
   # de eec
   nbcho[eec$IDENT[1:10]]
   ## G0A56JP6 G0A56JP6 G0A56JR6 G0A56JS6 G0A56JT6 G0A56JU6 G0A56JV6 G0A56JV6 
   ##        1        1        0        0        0        0        0        0 
   ## G0A56JW6 G0A56JX6 
   ##        0        0
   
   # On peut donc par ce biais reconstituer l'intégralité 
   # du vecteur correspondant aux observations de eec.
   eec$nbcho2 <- nbcho[eec$IDENT]
   
   # Ce vecteur est bien identique à celui obtenu par fusion
   # à la question précédente.
   all.equal(eec$nbcho1, eec$nbcho2)
   ## [1] TRUE

   ---

    

4. Comparez la syntaxe et les performances des méthodes mises en oeuvre aux question b. et c.

   Afficher/masquer la solution

   ---

   # La syntaxe de la deuxième option est plus concise
   # mais aussi plus complexe pour un relecteur moins
   # averti des fonctionnalités de R en matière d'utilisation
   # des noms de vecteur. 
   
   # Comparaison des performances
   microbenchmark(times = 10
     , merge1 = merge(eec, data.frame(IDENT = names(nbcho), nbcho1 = nbcho), by = "IDENT")
     , merge2 = merge(eec, nbcho_df, by = "IDENT")
     , names = nbcho[eec$IDENT]
   )
   ## Unit: milliseconds
   ##    expr        min         lq       mean     median         uq        max
   ##  merge1 111.404533 111.782110 119.345895 122.373338 124.426139 128.791746
   ##  merge2  94.860742  94.923989 100.040972  97.359094 105.536063 110.802111
   ##   names   2.284167   2.405333   2.492786   2.429073   2.472297   3.109997
   ##  neval
   ##     10
   ##     10
   ##     10
   # C'est sans commune mesure : que l'on tienne 
   # compte (merge1) ou pas (merge2) de l'étape
   # de constitution du data.frame, la méthode
   # reposant sur les noms de vecteur est beaucoup
   # plus rapide.

   ---

    

Travailler efficacement sur des données avec dplyr

Les cas pratiques de cette partie reposent sur le package dplyr :

install.packages("dplyr")
library(dplyr)

Plusieurs vignettes sont disponibles sur la page de documentation du package. Rstudio a également conçu un aide-mémoire (cheatsheet) traduit en français. N’hésitez pas à vous référer à ces documents pour répondre aux cas pratiques de cette partie.

 

Cas pratique 11 Sélection d’observations, de variables et tris

1. Utilisez le verbe filter() pour afficher les observations des femmes (SEXE == "2") actives occupées (ACTEU == "1") en Île-de-France (REG == "11"). Utilisez la syntaxe classique puis celle faisant appel à l’opérateur pipe %>%.

   Afficher/masquer la solution

   ---

   # La fonction filter() permet de ne pas avoir à répéter
   # le nom de la table et de pouvoir séparer les 
   # différentes clauses par des ,
   filter(eec, SEXE == "2", ACTEU == "1", REG == "11")
   eec %>% filter(SEXE == "2", ACTEU == "1", REG == "11")

   ---

    

2. Utilisez le verbe select() pour supprimer toutes les variables créées à la sous-partie précédente (trage1, trage2, trage3, trage4, nbcho1, nbcho2). Pensez à consulter les exemples de l’aide de select() pour l’utiliser au mieux.

   Afficher/masquer la solution

   ---

   # select() dispose de nombreuses fonctions outils
   # pour simplement sélectionner (ou exclure avec -)
   # des variables selon les caractères qu'elles contiennent.
   eec %>% 
     select(-starts_with("trage"), -contains("nbcho")) -> 
     eec
   # Note : utilisé à l'issue de plusieurs instruction, 
   # l'opérateur `->` permet d'assigner des valeurs 
   # à un objet situé à sa droite (et non à sa gauche comme `<-`)

   ---

    

3. Utilisez le verbe arrange() pour trier la table par région et identifiant de ménage croissants puis par numéro d’ordre dans le ménage décroissants.

   Afficher/masquer la solution

   ---

   # Il suffit d'indiquer la liste des variables sur
   # lesquelles trier, éventuellement avec la fonction
   # desc() quand l'ordre est décroissant.
   eec %>% 
     arrange(REG, IDENT, desc(NOI)) -> 
     eec

   ---

    

 

Cas pratique 12 Création de variables

1. Utilisez le verbe mutate() pour effectuer le recodage en classes d’âge présenté lors de la partie précédente.

   Afficher/masquer la solution

   ---

   # mutate() est analogue à la fonction de base R
   # transform() mais permet l'utilisation directe
   # des variables créées dans le même appel de fonction.
   eec %>% 
     mutate(
   age_num = as.numeric(AGE)
   , trage_dplyr = ifelse(age_num < 31, "15-30", ifelse(age_num < 61, "31-60", "61 et +"))
     ) ->
     eec

   ---

    

2. Comparez l’ergonomie et les performances de mutate() avec la méthode la plus efficace de base R.

   Afficher/masquer la solution

   ---

   # mutate() est relativement ergonomique dans la mesure
   # où elle permet de ne pas répéter le nom de la table,
   # de créer simultanément plusieurs variables et de 
   # réutiliser immédiatement les variables créées.
   
   # Comparaison des performances
   microbenchmark(times = 100
     , base = {
   t <- as.numeric(eec$AGE)
   eec$trage4 <- "15-30"
   eec$trage4[t > 30 & t < 61] <- "31-60"
   eec$trage4[t > 60] <- "61 et +"
     }
     , dplyr = {
   eec %>% 
     mutate(
       age_num = as.numeric(AGE)
       , trage_dplyr = ifelse(age_num < 31, "15-30", ifelse(age_num < 61, "31-60", "61 et +"))
     ) ->
     eec
     }
   )
   ## Unit: milliseconds
   ##   expr       min       lq      mean    median        uq       max neval
   ##   base  4.109121  4.14015  4.213202  4.157384  4.172201  9.602891   100
   ##  dplyr 24.094716 24.22197 25.254545 24.285097 24.406161 33.347734   100
   
   # Le gain en termes d'ergonomie de mutate() a donc 
   # un coût non-négligeable en termes de performances.

   ---

    

 

Cas pratique 13 Agrégation par groupes : salaire moyen par région

Comme dans le cas pratique correspondant de la partie précédente, l’objectif est d’estimer le salaire moyen par région, d’abord sans pondération puis pondéré par le poids de sondage EXTRI1613.

1. Utilisez la fonction summarise() pour calculer le salaire moyen non-pondéré et pondéré pour l’ensemble de la France métropolitaine. Intercalez ensuite la fonction group_by() pour ventiler ces calculs par région.

   Afficher/masquer la solution

   ---

   # On commence par calculer les moyennes
   # sur l'ensemble de la France métropolitaine
   # grâce à la fonction summarise().
   eec %>% 
     summarise(
   nonpond = mean(SALRED, na.rm = TRUE)
   , pond = weighted.mean(SALRED, EXTRI1613, na.rm = TRUE)
     )
   
   # Pour ventiler les calculs par région, il suffit
   # d'intercaler la fonction group_by()
   eec %>% group_by(REG) %>%
     summarise(
   nonpond = mean(SALRED, na.rm = TRUE)
   , pond = weighted.mean(SALRED, EXTRI1613, na.rm = TRUE)
     )

   ---

    

2. Comparez l’ergonomie et les performances des fonctions de dplyr avec la méthode la plus efficace de base R.

   Afficher/masquer la solution

   ---

   # La syntaxe est beaucoup plus ergonomique avec
   # les fonctions de dplyr, notamment du fait que
   # la ventilation par une ou plusieurs variables
   # ne nécessite qu'une adaptation minimale du code.
   
   # Comparaison des performances
   microbenchmark(times = 100
     , base = tapply(eec$SALRED, eec$REG, mean, na.rm = TRUE)
     , dplyr = eec %>% group_by(REG) %>% 
     summarise(SALRED = mean(SALRED, na.rm = TRUE))
   )
   ## Unit: milliseconds
   ##   expr      min       lq     mean   median       uq      max neval
   ##   base 2.106066 2.160475 2.363685 2.183966 2.219745 10.99048   100
   ##  dplyr 5.001960 5.102613 5.250645 5.144364 5.207435 13.67685   100
   # Ici dplyr est aussi rapide sinon plus que base R.
   # On gagne donc sur tous les tableaux.

   ---

    

 

Cas pratique 14 Fusion de tables : recodage de la PCS

La variable CSE de la table eec code la Profession et catégorie socio-professionnelle (PCS) au niveau 3 de la nomenclature Insee (cf. le site de l’Insee). On souhaite passer du niveau 3 au niveau 2. Le fichier pcs2003_c_n4_n1.dbf est la table de passage entre l’ensemble des niveaux de la nomenclature (de 1 à 4).

1. Utilisez le package foreign et la fonction read.dbf() pour lire le fichier pcs2003_c_n4_n1.dbf.

   Afficher/masquer la solution

   ---

   # La fonction read.dbf() du package foreign
   # permet de lire les fichier .dbf
   library(foreign)
   pcs <- read.dbf("pcs2003_c_n4_n1.dbf")

   ---

    

2. Utilisez le verbe distinct() pour restreindre la table aux observations distinctes pour les niveaux N2 et N3. Créez également la variable CSE, version caractère de N3.

   Afficher/masquer la solution

   ---

   # Ces deux opérations peuvent être effectuées
   # en une seule instruction avec des %>%
   pcs %>% 
     distinct(N2, N3) %>% 
     mutate(CSE = as.character(N3)) ->
     pcs

   ---

    

3. Utilisez le verbe left_join() pour fusionner la table eec avec la table de passage des PCS de niveau 3 à niveau 2.

   Afficher/masquer la solution

   ---

   # La syntaxe est rendue particulièrement claire 
   # (et proche de SQL) par l'utilisation des %>%
   eec %>% 
     left_join(pcs, by = "CSE") ->
     eec

   ---

    

4. Une solution purement vectorielle en base R n’aurait-elle pas également été possible ? Comparez les performances de dplyr avec cette solution alternative.

   Afficher/masquer la solution

   ---

   # On aurait aussi pu réutiliser le mécanisme
   # vu au cas pratique 13 en passant par un vecteur nommé
   pcs2 <- setNames(pcs$N2, pcs$N3)
   eec$N2_base <- pcs2[eec$CSE]
   all.equal(eec$N2, eec$N2_base)
   ## [1] TRUE
   
   # Comparaison des performances
   microbenchmark(times = 10
     , base = pcs2[eec$CSE]
     , dplyr = eec %>% left_join(pcs, by = "CSE")
   )
   ## Unit: milliseconds
   ##   expr       min        lq      mean    median       uq       max neval
   ##   base 819.91563 820.11250 820.33217 820.26625 820.4731 821.02741    10
   ##  dplyr  10.32031  10.37436  10.45911  10.44775  10.5158  10.71828    10
   # Il est possible que les perforamnces de la version 
   # en base R soit affectées par le grand nombre de NA.

   ---

    

Travailler efficacement sur des données avec data.table

Les cas pratiques de cette partie reposent sur le package data.table :

install.packages("data.table")
library(data.table)

On crée le data.table correspondant au data.frame eec :

eec_dt <- data.table(eec)

Plusieurs vignettes sont disponibles sur la page de documentation du package. N’hésitez pas à vous référer à ces documents pour répondre aux cas pratiques de cette partie.

 

Cas pratique 15 Sélection d’observations et tris

1. Utilisez l’argument i de [ pour afficher les observations des femmes (SEXE == "2") actives occupées (ACTEU == "1") en Île-de-France (REG == "11"). Quelle différence constatez-vous avec une sélection dans un data.frame ?

   Afficher/masquer la solution

   ---

   # Dans un data.frame et avec le [ de base R, 
   # il est nécessaire de répéter le nom de la table
   eec[eec$SEXE == "2" & eec$ACTEU == "1" & eec$REG == "11", ]
   
   # Ce n'est pas le cas dans un data.table
   eec_dt[SEXE == "2" & ACTEU == "1" & REG == "11", ]

   ---

    

2. Utilisez la fonction setkey() pour faire de SEXE, ACTEU et REG des clés pour eec_dt et utilisez-les pour reproduire la sélection de la question précédente. Comparez alors l’ergonomie et les performances d’une sélection d’observations :

   1. avec une clause logique dans un data.frame ;
   2. avec une clause logique en utilisant le verbe filter() de dplyr ;
   3. avec une clause logique dans un data.table ;
   4. avec un jeu de clés dans un data.table.
   Afficher/masquer la solution

   ---

   # La fonction setkey() permet de facilement créer
   # des clés pour un data.table donné.
   setkey(eec_dt, SEXE, ACTEU, REG)
   
   # La sélection sur la base de clés s'effetue
   # avec un argument sous la forme d'une liste
   # DANS L'ORDRE DES CLES
   eec_dt[list("2", "1", "11")]

   # Comparaison des performances
   microbenchmark(times = 100
     , base = eec[eec$SEXE == "2" & eec$ACTEU == "1" & eec$REG == "11", ]
     , dplyr = eec %>% filter(SEXE == "2", ACTEU == "1", REG == "11")
     , data.table1 = eec_dt[SEXE == "2" & ACTEU == "1" & REG == "11"]
     , data.table2 = eec_dt[list("2", "1", "11")]
   )
   ## Unit: milliseconds
   ##         expr      min       lq     mean   median       uq       max neval
   ##         base 3.674216 3.748339 4.377710 3.762984 3.792696 11.601206   100
   ##        dplyr 7.909269 8.000830 8.125151 8.043670 8.132972  9.688663   100
   ##  data.table1 1.575499 1.668389 1.723698 1.713931 1.761211  2.104813   100
   ##  data.table2 1.334407 1.421158 1.489781 1.504225 1.550642  2.036332   100
   
   # Les différentes version de data.table sont plus efficaces
   # que base R et dplyr, en particulier quand il est fait
   # usage des clés. 

   ---

    

3. Utilisez la fonction order() (comme dans un data.frame) pour trier la table eec_dt par région et identifiant de ménage croissants puis par numéro d’ordre dans le ménage décroissants. Comparez les performances de base R, arrange() de dplyr et data.table.

   Afficher/masquer la solution

   ---

   # La principale différence avec la fonction order()
   # appliquée à un data.frame est qu'il n'est pas 
   # nécessaire de répéter le nom de la table et
   # qu'il est possible d'utiliser le signe - devant 
   # des variables caractère
   eec_dt <- eec_dt[order(REG, IDENT, -NOI)]
   
   # Comparaison des performances
   microbenchmark(times = 10
     , base = eec[order(eec$REG, eec$IDENT, -as.numeric(eec$NOI)), ]
     , dplyr = eec %>% arrange(REG, IDENT, desc(NOI))
     , data.table = eec_dt[order(REG, IDENT, -NOI)]
   )
   ## Unit: milliseconds
   ##        expr      min       lq     mean   median       uq       max neval
   ##        base 90.19476 90.60685 92.74277 90.78227 91.04655 102.64125    10
   ##       dplyr 62.74580 62.97458 63.79840 63.78389 64.04781  65.20523    10
   ##  data.table 13.26530 13.38403 14.40743 13.42128 13.51141  23.14988    10
   # Le gain en termes de performances est sensible à nouveau.

   ---

    

 

Cas pratique 16 Agrégation par groupes : salaire moyen par région

Comme dans le cas pratiques correspondants des parties précédentes, on cherche à calculer le salaire moyen par région non-pondéré puis pondéré par le poids de sondage EXTRI1613.

1. Utilisez les arguments j et by de [ pour calculer le salaire non-pondéré et pondéré d’abord sur l’ensemble de la France métropolitaine, puis par région. Comparez l’utilisation de by et keyby : pourquoi les résultats sont-ils ici identiques à votre avis ?

   Afficher/masquer la solution

   ---

   # L'argument `j` de [ permet de créer de nouvelles 
   # variables dans un data.table
   eec_dt[, j = list(
     nonpond = mean(SALRED, na.rm = TRUE)
     , pond = weighted.mean(SALRED, EXTRI1613, na.rm = TRUE)
   )]
   ##     nonpond     pond
   ## 1: 1819.209 1833.879
   # Pour ventiler par région, il suffit d'ajouter un argument by
   eec_dt[, j = list(
     nonpond = mean(SALRED, na.rm = TRUE)
     , pond = weighted.mean(SALRED, EXTRI1613, na.rm = TRUE)
   ), by = REG]  
   ##     REG  nonpond     pond
   ##  1:  11 2153.045 2172.578
   ##  2:  21 1692.069 1686.691
   ##  3:  22 1585.062 1594.508
   ##  4:  23 1722.671 1727.632
   ##  5:  24 1702.413 1725.168
   ##  6:  25 1633.108 1625.475
   ##  7:  26 1640.106 1617.565
   ##  8:  31 1730.126 1714.465
   ##  9:  41 1729.073 1693.117
   ## 10:  42 1974.987 2012.017
   ## 11:  43 1747.486 1750.870
   ## 12:  52 1609.202 1616.999
   ## 13:  53 1784.060 1814.576
   ## 14:  54 1587.502 1547.045
   ## 15:  72 1703.824 1667.706
   ## 16:  73 1756.656 1829.931
   ## 17:  74 1633.896 1657.696
   ## 18:  82 1948.086 1954.849
   ## 19:  83 1633.669 1621.692
   ## 20:  91 1557.426 1538.386
   ## 21:  93 1804.827 1806.120
   ## 22:  94 1828.872 1796.574
   ##     REG  nonpond     pond
   # Les résultats sont identiques ici selon que l'on utilise
   # by ou keyby car les données sont triées par région.
   # Si cela n'était pas le cas, by conserverait l'ordre du
   # fichier (en affichant les groupes par ordre de rencontre)
   # alors que keyby trierait les résultats par ordre
   # alphabétique de région.

   ---

    

2. Comparez l’ergonomie et les performances de la solution en base R, avec dplyr et data.table.

   Afficher/masquer la solution

   ---

   # La solution en data.table est plus ergonomique que
   # celle en base R, en particulier en raison de la
   # facilité à ventiler les traitements par région. 
   # Néanmoins, elle ne dispose pas de la capacité 
   # à séquencer les traitements en petites opérations
   # simples, ce que permet l'opérateur %>% qu'utilise
   # intensément dplyr.
   
   # Comparaison des performances
   microbenchmark(times = 100
     , base = tapply(eec$SALRED, eec$REG, mean, na.rm = TRUE)
     , dplyr = eec %>% group_by(REG) %>% summarise(SALRED = mean(SALRED, na.rm = TRUE))
     , data.table = eec_dt[, j = list(mean(SALRED, na.rm = TRUE)), by = REG]
   )
   ## Unit: microseconds
   ##        expr      min       lq     mean   median       uq       max neval
   ##        base 2070.475 2151.584 2188.901 2172.514 2214.896  2589.206   100
   ##       dplyr 5066.877 5152.562 5589.764 5191.659 5241.400 14508.959   100
   ##  data.table  877.443  950.726 1159.415 1094.593 1117.881  9729.256   100
   # Là encore data.table est beaucoup plus rapide.

   ---

    

 

Cas pratique 17 Fusion de tables : nombre d’individus au chômage par ménage

Comme dans le cas pratique 13, on cherche ici à associer à chaque individu le nombre d’individus au chômage (ACTEU == "2") dans son ménage (individus avec la même valeur pour la variable IDENT).

1. Utilisez les arguments j et by de [ pour calculer le nombre d’individus au chômage par ménage. Utilisez la structure j := pour automatiquement refusionner ce résultat avec la table de départ.

   Afficher/masquer la solution

   ---

   # On reprend la syntaxe de la question précédente
   # pour calculer le nombre d'individu au chômage par ménage
   eec_dt[, sum(ACTEU == "2", na.rm = TRUE), by = "IDENT"]
   ##           IDENT V1
   ##     1: G0A56JP6  1
   ##     2: G0A56JR6  0
   ##     3: G0A56JS6  0
   ##     4: G0A56JT6  0
   ##     5: G0A56JU6  0
   ##    ---            
   ## 18899: GXZ5OX1F  0
   ## 18900: GY05O5DF  0
   ## 18901: GY05O85F  0
   ## 18902: GY05OAXF  1
   ## 18903: GY05ODPF  0
   
   # Pour refusionner ces résultats avec la table d'origine,
   # il suffit d'utiliser l'opérateur := au niveau de l'argument j
   eec_dt <- eec_dt[, nbcho := sum(ACTEU == "2", na.rm = TRUE), by = "IDENT"]

   ---

    

2. Comment mèneriez vous ce traitement dans la logique de dplyr ? Comparez l’ergonomie et les performances de la solution en base R, dplyr et data.table.

   Afficher/masquer la solution

   ---

   # Avec dplyr, on serait tenté de fusionner la table eec
   # avec une table de statistique (comme en base R)
   eec %>% 
     left_join(
   eec %>% group_by(IDENT) %>% summarize(nbcho = n())
   , by = "IDENT"
     ) -> 
     eec
   
   # Comparaison des performances
   microbenchmark(times = 10
     , base = tapply(eec$ACTEU == "2", eec$IDENT, sum, na.rm = TRUE)[eec$IDENT]
     , dplyr = {
   eec %>% 
     left_join(
       eec %>% group_by(IDENT) %>% summarize(nbcho = n())
       , by = "IDENT"
     )
     }
     , data.table = eec_dt[, nbcho := sum(ACTEU == "2", na.rm = TRUE), by = "IDENT"]
   )
   ## Unit: milliseconds
   ##        expr       min        lq      mean    median        uq       max
   ##        base  36.51546  36.59215  38.50903  36.71289  36.84654  45.88869
   ##       dplyr 123.16010 123.63943 125.27296 123.91513 124.27834 131.63315
   ##  data.table  16.24202  16.44046  17.46490  16.55652  16.71967  25.54495
   ##  neval
   ##     10
   ##     10
   ##     10
   # Comme toujours data.table est le plus rapide. 
   # Néanmoins ici, il n'est pas à exclure qu'il existe
   # dans dplyr une méthode plus efficace (une possibilité
   # d'autofusion après agrégation par exemple). 

   ---

    

Réaliser des graphiques avec R

Les cas pratiques de cette partie reposent sur l’utilisation du package ggplot2 :

install.packages("ggplot2")
library(ggplot2)

 

Cas pratique 18 Graphiques à partir de la table mpg

L’objectif de ce cas pratique est de reproduire les graphiques du support ainsi que ceux présentés par H. Wickham dans le chapitre 2 de son ouvrage ggplot2: Elegant Graphics for Data Analysis.

Dans les deux cas, la table utilisée est mpg (table d’exemple du package ggplot2) : après avoir chargé ggplot2, utilisez la fonction data() pour “rapatrier” la table mpg dans l’environnement global et tapez ? mpg pour obtenir une description détaillée de ses variables. En particulier :

1. Utilisez les mots-clés colour, shape et size pour faire varier la représentation des points avec geom_point(). Pour chacun des mots-clés, comparez ce qu’il se passe quand vous utilisez une variable de type numérique ou une variable de type caractère ou facteur.
2. Comparez l’utilisation du mot-clé colour dans la fonction aes() et en dehors de la fonction aes().
3. Testez les différents types de représentation possibles en les adaptant à la nature des données à représenter. Pour chaque fonction geom_*(), recherchez dans l’aide les paramètres qui lui sont spécifiques et testez des valeurs différentes de leurs valeurs par défaut.
4. Expérimentez les différentes possibilités de facetting.
5. Sauvegarder un graphique dans un objet R. Utilisez ggsave() pour exporter un graphique en .png et .pdf.
6. Affichez le code d’une fonction geom_*() et utilisez ces informations pour reconstituer manuellement l’instruction layer() correspondante.
7. Tentez de reproduire un graphique standard de ggplot2 avec les fonctions du package graphics.

 

Cas pratique 19 Graphiques à partir de la table diamonds

La table diamonds est le second fichier de démonstration classique de ggplot2 : utilisez data() pour le “rapatrier” dans l’environnement global et tapez ? diamonds pour obtenir une description détaillée de ses variables.

1. Représentez la relation entre poids du diamant (carat) et prix (price). Utilisez le paramètre alpha pour limiter la saturation du graphique par le très grand nombre de points. Ajoutez une droite de régression linéaire au graphique.

   Afficher/masquer la solution

   ---

   # Par défaut, les points sont totalement opaques : 
   # on ne peut pas visualiser les points qui se superposent
   # les uns aux autres (on parle d'over-plotting)
   ggplot(diamonds, aes(carat, price)) + geom_point()

   

   
   # Le paramètre alpha indique de rendre les points
   # en partie transparent. Avec alpha = 0.05, il faut
   # 20 points pour obtenir une zone totalement opaque
   ggplot(diamonds, aes(carat, price)) + geom_point(alpha = 0.05)

   

   
   # Pour ajouter une droite de régression, il suffit 
   # d'utiliser la fonction geom_smooth() avec l'option
   # method = "lm"
   ggplot(diamonds, aes(carat, price)) + geom_point(alpha = 0.05) + 
     geom_smooth(method = "lm")

   

   ---

    

2. Représentez l’influence de la couleur (color) sur le prix de plusieurs manières.

   Afficher/masquer la solution

   ---

   # Idée 1 : faire varier la couleur des points 
   # sur le graph précédent
   g1 <- ggplot(diamonds, aes(carat, price, colour = color)) + geom_point(alpha = 0.05)
   g1

   

   
   # Il est possible d'améliorer le graphique précédent
   # 1) en ajustant la manière dont les couleurs sont représentées 
   # dans la légende
   g1 <- g1 + guides(colour = guide_legend(override.aes = list(alpha = 1)))
   g1

   

   
   # 2) en adoptant une palette de couleurs formant un gradient
   # (car la variable color est ordonnée de la pire (J) à la 
   # meilleure (D))
   g1 <- g1 + scale_colour_brewer(palette = 1)
   g1

   

   
   # Idée 2 : dessiner des boîtes à moustaches
   g2 <- ggplot(diamonds, aes(color, price)) + geom_boxplot(varwidth = TRUE)
   g2

   

   ---

    

3. Représentez la ventilation des diamants selon la qualité de leur taille (cut) et leur clarté (clarity).

   Afficher/masquer la solution

   ---

   # Le plus simple est de représenter l'histogramme bivarié
   ggplot(diamonds, aes(clarity, fill = cut)) + geom_bar()

   

   
   # Quelques variations sur le positionnement des blocs
   ggplot(diamonds, aes(clarity, fill = cut)) + geom_bar(position = "dodge")

   

   ggplot(diamonds, aes(clarity, fill = cut)) + geom_bar(position = "fill")

   

   ---

    

4. Vérifiez graphiquement si la relation entre poids et prix ne varie pas en fonction de la qualité de la taille (cut) et la clarté du diamant (clarity).

   Afficher/masquer la solution

   ---

   # L'idée ici est d'utiliser le facetting pour ventiler
   # le premier graphique par qualité de la taille et clarté 
   g3 <- ggplot(diamonds, aes(carat, price)) + geom_smooth() +
     facet_grid(cut ~ clarity)
   g3
   ## `geom_smooth()` using method = 'gam'
   ## Warning: Computation failed in `stat_smooth()`:
   ## x has insufficient unique values to support 10 knots: reduce k.

   

   ---

    

 

Cas pratique 20 Graphiques à partir de la table raisin

Le fichier raisin.rds comporte des informations sur la maturation du raisin dans des exploitations viticoles de Saône-et-Loire sur la période 2000-2012.

1. Chargez ce fichier en mémoire avec la fonction readRDS() et analysez les caractéristiques des variables de ce fichier (modalités, distributions, etc.).

   Afficher/masquer la solution

   ---

   # Chargement en mémoire du fichier raisin.rds
   # situé dans le répertoire de travail
   raisin <- readRDS("raisin.rds")
   
   # Caractéristiques des variables de raisin
   table(raisin$annee)
   ## 
   ## 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 
   ##   54   48   56   39   76   72   79  124  116  111  124  117  111
   table(raisin$secteur)
   ## 
   ##  Couchois Maconnais  Maranges 
   ##       110       907       110
   table(raisin$cepage)
   ## 
   ## Chardonnay      Gamay      Pinot 
   ##        479        428        220
   table(raisin$commune)
   ## 
   ##   Cheilly-les-Maranges                 Fuissé               St.Amour 
   ##                    110                    479                    428 
   ## St.Maurice-les-Couches     St.Sernin-du-Plain 
   ##                     55                     55
   summary(raisin)
   ##    secteur             cepage            commune              sucres      
   ##  Length:1127        Length:1127        Length:1127        Min.   :  48.0  
   ##  Class :character   Class :character   Class :character   1st Qu.: 136.0  
   ##  Mode  :character   Mode  :character   Mode  :character   Median : 167.0  
   ##                                                           Mean   : 348.4  
   ##                                                           3rd Qu.: 197.0  
   ##                                                           Max.   :2057.0  
   ##                                                                           
   ##  acidite_totale         ph        acide_tartrique    potasse     
   ##  Min.   : 3.000   Min.   : 1.00   Min.   : 1.00   Min.   : 1.00  
   ##  1st Qu.: 7.000   1st Qu.: 4.00   1st Qu.: 9.00   1st Qu.:12.00  
   ##  Median : 8.000   Median : 7.00   Median :74.00   Median :19.00  
   ##  Mean   : 8.987   Mean   :21.56   Mean   :55.91   Mean   :21.22  
   ##  3rd Qu.:11.000   3rd Qu.:32.00   3rd Qu.:87.00   3rd Qu.:26.00  
   ##  Max.   :21.000   Max.   :99.00   Max.   :99.00   Max.   :99.00  
   ##                                                   NA's   :1      
   ##    azote_amm         annee     
   ##  Min.   : 1.00   Min.   :2000  
   ##  1st Qu.:12.00   1st Qu.:2005  
   ##  Median :43.00   Median :2008  
   ##  Mean   :43.76   Mean   :2007  
   ##  3rd Qu.:71.00   3rd Qu.:2010  
   ##  Max.   :99.00   Max.   :2012  
   ##  NA's   :28

   ---

    

2. Analysez la fréquence des différents cépages en fonction du temps.

   Afficher/masquer la solution

   ---

   # Le plus simple est ici de faire un diagramme en bâton
   # en fonction du temps
   ggplot(raisin, aes(annee, fill = cepage)) + geom_bar()

   

   
   # On peut utiliser explicitement la statistique "count"
   # associée par défaut à geom_bar() (tapez geom_bar pour
   # le vérifier) à d'autres fonctions pour modifier cette
   # représentation
   ggplot(raisin, aes(annee, colour = cepage)) + geom_line(stat = "count")

   

   ggplot(raisin, aes(annee, fill = cepage)) + geom_area(stat = "count")

   

   ---

    

3. Analysez graphiquement la relation entre sucres et acidite_totale. Utilisez d’autres variables pour tenter de rendre compte de cette distribution.

   Afficher/masquer la solution

   ---

   # On utilise tout simplement un nuage de points pour représenter
   # ces deux variables quantitatives
   ggplot(raisin, aes(sucres, acidite_totale)) + geom_point()

   

   
   # On perçoit très nettement deux groupes : sont-ils expliqués
   # par le secteur ou le cépage ? 
   ggplot(raisin, aes(sucres, acidite_totale, colour = cepage)) + geom_point()

   

   ggplot(raisin, aes(sucres, acidite_totale, colour = secteur)) + geom_point()

   

   
   # Pour confirmer et ne pas se laisser abuser
   # par la superposition de certains points, on 
   # utilise le facetting
   ggplot(raisin, aes(sucres, acidite_totale)) + geom_point() + 
     facet_wrap(~cepage)

   

   ggplot(raisin, aes(sucres, acidite_totale)) + geom_point() + 
     facet_wrap(~secteur)

   

   ggplot(raisin, aes(sucres, acidite_totale)) + geom_point() + 
     facet_grid(cepage~secteur)

   

   # Il apparaît clairement que c'est la localisation dans le
   # maconnais ou les cépages chardonnay et gamay qui lui
   # sont spécifiques (au sein de la Saône-et-Loire) qui 
   # semblent expliquer l'oscillation entre deux types de relations
   # entre sucres et acidité totale. 
   
   # Peut-être cette oscillation dépend-elle des années ? 
   ggplot(raisin[raisin$secteur == "Maconnais", ], aes(sucres, acidite_totale, colour = as.factor(annee))) + geom_point()

   

   
   # Ce n'est pas évident, à nouveau on utilise le facetting
   ggplot(raisin[raisin$secteur == "Maconnais", ], aes(sucres, acidite_totale, colour = commune)) + geom_point() + facet_wrap(~annee)

   

   # Le raisin du secteur du maconnais, quel que soit sa commune, 
   # semble présenter un profil sucres-acidite_totale particulier
   # de 2000 à 2003 (plus sucré).

   ---