Ce chapitre vise à illustrer l’utilisation de la notation formule de R, qui désigne l’emploi de cette notation par l’expression formula. Cette notation est utilisée par de très nombreuses fonctions de R : on en a notamment vu plusieurs exemples dans le chapitre sur les graphiques bivariés, car l’extension ggplot2 se sert de cette notation dans ses paramètres facet_wrap et facet_grid.

Dans ce chapitre, on verra comment se servir de la notation formule dans deux contextes différents. D’une part, on verra que deux fonctions basiques de R se servent de cette notation pour produire des tableaux croisés et des statistiques bivariées. D’autre part, on verra que l’extension lattice se sert de cette notation pour créer des graphiques panelisés, dits graphiques à petits multiples.

Dans plusieurs autres chapitres, les opérations décrites ci-dessus sont effectuées avec les extensions dplyr d’une part, et ggplot2 d’autre part. On se servira également de ces extensions dans ce chapitre, de manière à mener une comparaison des différentes manières d’effectuer certaines opérations dans R, avec ou sans la notation formule :

library(dplyr)
library(ggplot2)

Statistiques descriptives

Les premiers exemples de ce chapitre montrent l’utilisation de cette notation pour produire des tableaux croisés et des statistiques descriptives. Le jeu de données utilisé, hdv2003, a déjà été utilisé dans plusieurs chapitres, et font partie de l’extension questionr. Chargeons cette extension et le jeu de données hdv2003 :

library(questionr)
data(hdv2003)

Pour rappel, ce jeu de données contient des individus, leur âge, leur statut professionnel, et le nombre d’heures quotidiennes passées à regarder la télévision.

glimpse(hdv2003, 75)
Observations: 2,000
Variables: 20
$ id            <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 1...
$ age           <int> 28, 23, 59, 34, 71, 35, 60, 47, 20, 28, 65, 47, ...
$ sexe          <fctr> Femme, Femme, Homme, Homme, Femme, Femme, Femme...
$ nivetud       <fctr> Enseignement superieur y compris technique supe...
$ poids         <dbl> 2634.3982, 9738.3958, 3994.1025, 5731.6615, 4329...
$ occup         <fctr> Exerce une profession, Etudiant, eleve, Exerce ...
$ qualif        <fctr> Employe, NA, Technicien, Technicien, Employe, E...
$ freres.soeurs <int> 8, 2, 2, 1, 0, 5, 1, 5, 4, 2, 3, 4, 1, 5, 2, 3, ...
$ clso          <fctr> Oui, Oui, Non, Non, Oui, Non, Oui, Non, Oui, No...
$ relig         <fctr> Ni croyance ni appartenance, Ni croyance ni app...
$ trav.imp      <fctr> Peu important, NA, Aussi important que le reste...
$ trav.satisf   <fctr> Insatisfaction, NA, Equilibre, Satisfaction, NA...
$ hard.rock     <fctr> Non, Non, Non, Non, Non, Non, Non, Non, Non, No...
$ lecture.bd    <fctr> Non, Non, Non, Non, Non, Non, Non, Non, Non, No...
$ peche.chasse  <fctr> Non, Non, Non, Non, Non, Non, Oui, Oui, Non, No...
$ cuisine       <fctr> Oui, Non, Non, Oui, Non, Non, Oui, Oui, Non, No...
$ bricol        <fctr> Non, Non, Non, Oui, Non, Non, Non, Oui, Non, No...
$ cinema        <fctr> Non, Oui, Non, Oui, Non, Oui, Non, Non, Oui, Ou...
$ sport         <fctr> Non, Oui, Oui, Oui, Non, Oui, Non, Non, Non, Ou...
$ heures.tv     <dbl> 0.0, 1.0, 0.0, 2.0, 3.0, 2.0, 2.9, 1.0, 2.0, 2.0...

Tableaux croisés avec xtabs

Utilisons, pour ce premier exemple, la variable occup du jeu de données hdv2003, qui correspond au statut professionnel des individus inclus dans l’échantillon. La fonction de base pour compter les individus par statut est la fonction table :

table(hdv2003$occup)

Exerce une profession               Chomeur 
                 1049                   134 
      Etudiant, eleve              Retraite 
                   94                   392 
  Retire des affaires              Au foyer 
                   77                   171 
        Autre inactif 
                   83 

Avec la fonction xtabs, le même résultat est produit à partir de la notation suivante :

xtabs(~occup, data = hdv2003)
occup
Exerce une profession               Chomeur 
                 1049                   134 
      Etudiant, eleve              Retraite 
                   94                   392 
  Retire des affaires              Au foyer 
                   77                   171 
        Autre inactif 
                   83 

Le premier argument est une formule, au sens où R entend cette expression. Le second argument, data, correspond au jeu de données auquel la formule doit être appliquée. On pourra se passer d’écrire explicitement cet argument dans les exemples suivants.

L’avantage de la fonction xtabs n’est pas évident dans ce premier exemple. En réalité, cette fonction devient utile lorsque l’on souhaite construire un ou plusieurs tableau(x) croisé(s). Par exemple, pour croiser la variable occup avec la variable sexe, une solution constiste à écrire :

with(hdv2003, table(occup, sexe))
                       sexe
occup                   Homme Femme
  Exerce une profession   520   529
  Chomeur                  54    80
  Etudiant, eleve          48    46
  Retraite                208   184
  Retire des affaires      39    38
  Au foyer                  0   171
  Autre inactif            30    53

Ou alors, ce qui revient au même :

table(hdv2003$occup, hdv2003$sexe)

Avec xtabs, la même opération s’écrit de la manière suivante :

xtabs(~occup + sexe, hdv2003)
                       sexe
occup                   Homme Femme
  Exerce une profession   520   529
  Chomeur                  54    80
  Etudiant, eleve          48    46
  Retraite                208   184
  Retire des affaires      39    38
  Au foyer                  0   171
  Autre inactif            30    53

Cette écriture est plus courte que le code équivalent dans dplyr :

group_by(hdv2003, occup) %>%
  summarise(Homme = sum(sexe == "Homme"),
            Femme = sum(sexe == "Femme"))
# A tibble: 7 × 3
                  occup Homme Femme
                 <fctr> <int> <int>
1 Exerce une profession   520   529
2               Chomeur    54    80
3       Etudiant, eleve    48    46
4              Retraite   208   184
5   Retire des affaires    39    38
6              Au foyer     0   171
7         Autre inactif    30    53

De plus, xtabs permet de créer plusieurs tableaux croisés en une seule formule :

xtabs(~occup + sexe + trav.imp, hdv2003)
, , trav.imp = Le plus important

                       sexe
occup                   Homme Femme
  Exerce une profession    13    16
  Chomeur                   0     0
  Etudiant, eleve           0     0
  Retraite                  0     0
  Retire des affaires       0     0
  Au foyer                  0     0
  Autre inactif             0     0

, , trav.imp = Aussi important que le reste

                       sexe
occup                   Homme Femme
  Exerce une profession   159   100
  Chomeur                   0     0
  Etudiant, eleve           0     0
  Retraite                  0     0
  Retire des affaires       0     0
  Au foyer                  0     0
  Autre inactif             0     0

, , trav.imp = Moins important que le reste

                       sexe
occup                   Homme Femme
  Exerce une profession   328   380
  Chomeur                   0     0
  Etudiant, eleve           0     0
  Retraite                  0     0
  Retire des affaires       0     0
  Au foyer                  0     0
  Autre inactif             0     0

, , trav.imp = Peu important

                       sexe
occup                   Homme Femme
  Exerce une profession    20    32
  Chomeur                   0     0
  Etudiant, eleve           0     0
  Retraite                  0     0
  Retire des affaires       0     0
  Au foyer                  0     0
  Autre inactif             0     0

Cet exemple permet simplement de réaliser que la variable trav.imp, qui contient les réponses à une question portant sur l’importance du travail, n’a été mesurée (c’est-à-dire que la question n’a été posée) qu’aux seuls individus actifs de l’échantillon.

Statistiques bivariées avec aggregate

aggregate(heures.tv ~ sexe, mean, data = hdv2003)
   sexe heures.tv
1 Homme  2.219330
2 Femme  2.268727

Ici, le premier argument est à nouveau une formule. Le second argument correspond à la statistique descriptive que l’on souhaite obtenir, et le dernier argument indique le jeu de données auquel appliquer les deux autres arguments. On peut d’ailleurs obtenir le même résultat en respectant de manière plus stricte l’ordre des arguments dans la syntaxe de la fonction aggregate :

aggregate(heures.tv ~ sexe, hdv2003, mean)
   sexe heures.tv
1 Homme  2.219330
2 Femme  2.268727

Cette écriture est, à nouveau, plus compacte que le code équivalent dans dplyr, qui demande de spécifier le retrait des valeurs manquantes :

group_by(hdv2003, sexe) %>%
  summarise(heures.tv = mean(heures.tv, na.rm = TRUE))

À nouveau, on va pouvoir combiner plusieurs variables dans la formule que l’on passe à aggregate, ce qui va permettre d’obtenir la moyenne des heures de télévision quotidiennes par sexe et par statut professionnel :

aggregate(heures.tv ~ sexe + occup, hdv2003, mean)
    sexe                 occup heures.tv
1  Homme Exerce une profession  1.920463
2  Femme Exerce une profession  1.724953
3  Homme               Chomeur  2.853846
4  Femme               Chomeur  2.888608
5  Homme       Etudiant, eleve  1.400000
6  Femme       Etudiant, eleve  1.256522
7  Homme              Retraite  2.826442
8  Femme              Retraite  2.877174
9  Homme   Retire des affaires  2.410256
10 Femme   Retire des affaires  2.844737
11 Femme              Au foyer  2.822222
12 Homme         Autre inactif  3.133333
13 Femme         Autre inactif  3.339623

La même opération demanderait toujours un peu plus de code avec dplyr :

group_by(hdv2003, occup, sexe) %>%
  summarise(heures.tv = mean(heures.tv, na.rm = TRUE))

La fonction aggregate permet bien sûr d’utiliser une autre fonction que la moyenne, comme dans cet exemple, suivi de son équivalent avec dplyr :

# âge médian par sexe et statut professionnel
aggregate(age ~ sexe + occup, hdv2003, median)
    sexe                 occup  age
1  Homme Exerce une profession 42.0
2  Femme Exerce une profession 41.0
3  Homme               Chomeur 35.5
4  Femme               Chomeur 40.0
5  Homme       Etudiant, eleve 20.0
6  Femme       Etudiant, eleve 20.5
7  Homme              Retraite 68.0
8  Femme              Retraite 69.0
9  Homme   Retire des affaires 70.0
10 Femme   Retire des affaires 74.0
11 Femme              Au foyer 51.0
12 Homme         Autre inactif 56.0
13 Femme         Autre inactif 58.0
# code équivalent avec l'extension 'dplyr'
group_by(hdv2003, occup, sexe) %>%
  summarise(age = median(age, na.rm = TRUE))

La fonction aggregate permet, par ailleurs, d’obtenir des résultats à plusieurs colonnes. Dans l’exemple ci-dessus, on illustre ce principe avec la fonction range, qui renvoie deux résultats (la valeur minimale et la valeur maximale de la variable, qui est toujours la variable age), chacun présentés dans une colonne :

aggregate(age ~ sexe + occup, hdv2003, range)
    sexe                 occup age.1 age.2
1  Homme Exerce une profession    18    63
2  Femme Exerce une profession    18    67
3  Homme               Chomeur    18    63
4  Femme               Chomeur    18    63
5  Homme       Etudiant, eleve    18    34
6  Femme       Etudiant, eleve    18    35
7  Homme              Retraite    48    92
8  Femme              Retraite    41    96
9  Homme   Retire des affaires    57    91
10 Femme   Retire des affaires    57    93
11 Femme              Au foyer    22    90
12 Homme         Autre inactif    39    71
13 Femme         Autre inactif    19    97

Cette fonction ne peut pas être facilement écrite dans dplyr sans réécrire chacune des colonnes, ce que le bloc de code suivant illustre. On y gagne en lisibilité dans les intitulés de colonnes :

group_by(hdv2003, occup, sexe) %>%
  summarise(min = min(age, na.rm = TRUE),
            max = max(age, na.rm = TRUE))

Panels graphiques avec lattice

Les exemples suivants montreront ensuite comment la notation formule peut servir à produire des graphiques par panel avec l’extension lattice.

library(lattice)

L’extension lattice présente l’avantage d’être installée par défaut avec R. Il n’est donc pas nécessaire de l’installer préalablement.

Chargeons les mêmes données que le chapitre sur les graphiques bivariés.

# charger l'extension lisant le format CSV
library(readr)

# emplacement souhaité pour le jeu de données
file = "data/debt.csv"

# télécharger le jeu de données s'il n'existe pas
if(!file.exists(file))
  download.file("http://www.stat.cmu.edu/~cshalizi/uADA/13/hw/11/debt.csv",
                file, mode = "wb")

# charger les données dans l'objet 'debt'
debt = read_csv(file)
Warning: Missing column names filled in: 'X1' [1]
Parsed with column specification:
cols(
  X1 = col_integer(),
  Country = col_character(),
  Year = col_integer(),
  growth = col_double(),
  ratio = col_double()
)

Rejetons rapidement un coup d’oeil à ces données, qui sont structurées par pays (variable CountryYear). On y trouve deux variables, growth (le taux de croissance du produit intérieur brut réel), et ratio (le ratio entre la dette publique et le produit intérieur brut), ainsi qu’une première colonne vide, ne contenant que des numéros lignes, dont on va se débarrasser :

# inspection des données
glimpse(debt, 75)
Observations: 1,171
Variables: 5
$ X1      <int> 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157,...
$ Country <chr> "Australia", "Australia", "Australia", "Australia", "A...
$ Year    <int> 1946, 1947, 1948, 1949, 1950, 1951, 1952, 1953, 1954, ...
$ growth  <dbl> -3.5579515, 2.4594746, 6.4375341, 6.6119938, 6.9202012...
$ ratio   <dbl> 190.41908, 177.32137, 148.92981, 125.82870, 109.80940,...
# suppression de la première colonne
debt = debt[, -1]

Visualisation bivariée

Le même graphique s’écrit de la manière suivante avec l’extension lattice :

xyplot(growth ~ Year, data = debt)

Visualisation par petits multiples

Appliquons désormais la même visualisation par petits multiples que vue dans le chapitre :

xyplot(growth ~ Year | Country, data = debt)

Enfin, rajoutons quelques options au graphique, afin de montrer comment l’extension lattice fonctionne :

xyplot(growth ~ Year | Country, type = c("o", "l"), 
  main = "Données Reinhart et Rogoff corrigées, 1946-2009", 
  ylab = "Taux de croissance du PIB", xlab = NULL, 
  data = debt)

Pour aller plus loin

Comme vient de le voir dans ce chapitre, la notation formule apparaît çà et là dans les différentes fonctions de R est de ses extensions. Il est par conséquent utile d’en connaître les rudiments, et en particulier les opérateurs ~ (tilde) et +, ne serait-ce que pour pouvoir se servir des différentes fonctions présentées sur cette page.

La notation formule devient cruciale dès que l’on souhaite rédiger des modèles : la formule y ~ x, par exemple, qui est équivalente à la formule y ~ 1 + x, correspond à l’équation mathématique Y = a + bX. On trouvera de nombreux exemples d’usage de cette notation dans les chapitres consacrés, notamment, à la régression linéaire ou à la régression logistique.

De la même manière, l’opérateur | (pipe) utilisé par l’extension lattice joue aussi un rôle très important dans la rédaction de modèles multi-niveaux, où il sert à indiquer les variables à pentes ou à coefficients aléatoires. Ces modèles sont présentés dans un chapitre dédié.