Système graphique de base

Données quakes

Plusieurs des exercices suivants vous feront manipuler les données stockées dans le data frame quakes, inclus dans le package datasets de l’installation de base de R. Il s’agit d’observations pour 1000 évènements sismiques survenus près des îles Fidji depuis 1964. Ce sont les données utilisées dans les notes de cours.

Vous devez préliminairement ajouter à ces données deux facteurs nommés mag_catego et region en soumettant les commandes suivantes.

quakes$mag_catego <- factor(floor(quakes$mag))
quakes$region <- factor(
  ifelse(quakes$long < 175, yes = "Ouest", no = "Est"), 
  levels = c("Ouest", "Est")
)

Question 1 :

Fonction générique plot

Que produit la fonction plot lorsqu’elle reçoit en entrée comme premier argument

  1. un objet de classe ts, tel que l’objet suivant
WWWusage
  1. un objet de classe table à une ou deux dimensions, tel que les objets suivants
freq1 <- table(quakes$mag_catego)
freq2 <- table(quakes$mag_catego, quakes$region)

Solution :

plot(WWWusage)
# Diagramme en lignes ou temporel

freq1 <- table(quakes$mag_catego)
plot(freq1)
# Fréquences simples =>
# Genre de diagramme en bâtons pour lequel les bâtons ne sont que des lignes
# partant de 0 et montant verticalement jusqu'à la fréquence.

freq2 <- table(quakes$mag_catego, quakes$region)
plot(freq2)
# Fréquences croisées =>
# Diagramme en mosaïque

Question 2 :

Arguments graphiques avec le système graphique de base

Avec la fonction hist, créez un histogramme des observations de la variable stations du jeu de données quakes.

Assurez-vous d’avoir des annotations textuelles informatives et en français.

Limitez votre graphique aux observations entre 0 et 100 inclusivement.

Solution :

hist(quakes$stations)

# Comment modifier le titre ?
hist(
  quakes$stations,
  main = "Distribution du nombre de stations ayant rapporté un séisme"
)

# Comment modifier les noms des axes?
hist(
  quakes$stations,
  main = "Distribution du nombre de stations ayant rapporté un séisme",
  xlab = "nombre de stations", ylab = "fréquence"
)
# Comment modifier l'étendue de l'axe des x ?
# Nous voulons laisser tomber les valeurs au-delà de 100.
hist(
  quakes$stations,
  main = "Distribution du nombre de stations ayant rapporté un séisme",
  xlab = "nombre de stations", ylab = "fréquence",
  xlim = c(0, 100)
)
# Le résultat obtenu ne convient pas (dernier bâton partiellement visible).

# Ici, le résultat est mieux en donnant en entrée un sous-ensemble des données.
hist(
  quakes$stations[quakes$stations <= 100],
  main = "Distribution du nombre de stations ayant rapporté un séisme",
  xlab = "nombre de stations", ylab = "fréquence"
)

Question 3 :

Solutions pour observations superposées avec le système de base

Dans le diagramme de dispersion suivant

plot(
  x = quakes$mag, y = quakes$stations,
  xlab = "magnitude du séisme", 
  ylab = "nombre de stations rapportant le séisme"
)

la majorité des points représentent en fait plusieurs observations.

Expérimentez les deux solutions suivantes pour mieux visualiser les observations superposées. Ces solutions sont similaires à celles réalisées avec ggplot2 présentées dans les notes de cours à la section : https://stt4230.rbind.io/communication_resultats/graphiques_ggplot2_r/#solutions-aux-observations-superpos%C3%A9es-overplotting

  1. ajoutez un peu de bruit aux observations (sur les deux axes), à l’aide de la fonction jitter.

  2. produisez le diagramme de dispersion avec la fonction smoothScatter,

Solution :

# Afin de savoir comment utiliser la fonction jitter :
# help(jitter)

# Il faut transformer les observations avec jitter,
# puis les donner en entrée à la fonction plot.
plot(
  x = jitter(quakes$mag, factor = 2), 
  y = jitter(quakes$stations, factor = 2),
  xlab = "magnitude du séisme", 
  ylab = "nombre de stations rapportant le séisme"
)
# L'argument factor = 2 n'est pas obligatoire,
# mais je trouve que ça donne un meilleur résultat ici.
# Plus la valeur fournie à factor est élevée, plus la variabilité
# du bruit ajouté est grande.

# Afin de savoir comment utiliser la fonction smoothScatter : 
# help(smoothScatter)

# Seulement deux arguments obligatoires : x et y
smoothScatter(
  x = quakes$mag, y = quakes$stations,
  xlab = "magnitude du séisme", 
  ylab = "nombre de stations rapportant le séisme"
)

Question 4 :

Ajout d’éléments à un graphique produit avec le système graphique de base

Comment ajouter au diagramme de dispersion initial de la question précédente :

  1. la courbe de régression suivante
outreg <- lm(stations ~ mag, data = quakes)
  1. la courbe de lissage suivante
outlowess <- lowess(quakes$mag, quakes$stations)
  1. une légende, en haut à gauche, pour identifier les courbes ajoutées;

  2. l’erreur standard résiduelle (residual standard error) de la régression, en bas à droite, avec le libellé \(\sigma_e = 11.5\) ?

outregRSE <- round(summary(outreg)$sigma, 1)
outregRSE
## [1] 11.5

Solution :

plot(
  stations ~ mag, data = quakes, 
  xlab = "magnitude du séisme", 
  ylab = "nombre de stations rapportant le séisme",
  col = grey(0.5)
)

# i. Ajout de la courbe de régression
outreg <- lm(stations ~ mag, data = quakes)
abline(outreg, lty = 2)

# ii. Ajout de la courbe de lissage
outlowess <- lowess(quakes$mag, quakes$stations)
lines(outlowess, lty = 1)

# iii. Ajout de la légende
legend(
  x = "topleft", lty = c(2, 1), bty = "n",
  legend = c("droite de régression", "courbe de lissage")
)

# iv. Ajout textuel
outregRSE <- round(summary(outreg)$sigma, 1)
# text(x = 6, y = 15, label = expression(sigma[e] == 11.5))
# ou mieux
text(
  x = 6, y = 15, 
  label = substitute(
    exp = paste(sigma[e], " = ", valeur), 
    env = list(valeur = outregRSE)
  )
)

Question 5 :

Sous-fenêtres graphiques avec le système graphique de base

Reproduisez le graphique suivant, composé de plusieurs sous-graphiques, représentant des variables du jeu de données quakes.

 

Solution :

# Configuration de la disposition des sous-fenêtres graphiques
ng <- layout(
  matrix(
    c(2, 0, 3, 
      0, 1, 0, 
      4, 0, 5), 
    nrow = 3, byrow = TRUE
  ), 
  widths = c(1,2,1),
  heights = c(1,2,1)
)

# Pour visualiser la disposition des sous-fenêtres graphiques
# layout.show(n = ng)

# Ajustement des marges graphiques pour améliorer l'allure du résultat
par.default <- par(mar = c(4, 4, 0, 0) + 0.1)

# Graphique 1 (au centre)
plot(stations ~ mag, data = quakes)
# Graphique 2 (en haut à gauche)
boxplot(quakes$stations, ylab = "stations")
# Graphique 3 (en haut à droite)
boxplot(stations ~ mag_catego, data = quakes, xlab = "mag_catego", ylab = "stations")
# Graphique 4 (en bas à gauche)
boxplot(quakes$mag, xlab = "mag", horizontal = TRUE)
# Graphique 5 (en bas à droite)
barplot(table(quakes$mag_catego), xlab = "mag_catego", ylab = "Fréquence")

# Réattribution des valeurs par défaut aux paramètres graphiques
par(par.default)
layout(matrix(1))

Gaphiques avec ggplot2

Chargement du package ggplot2 :

library(ggplot2)

Question 6 :

Mapping et Facets avec ggplot2

Reproduisez le graphique suivant, composé de plusieurs sous-graphiques selon les niveaux des variables mag_catego et region, représentant des variables du jeu de données quakes.

Solution :

ggplot(data = quakes) +
  geom_point(
    mapping = aes(x = long, y = lat), 
    alpha = 0.3  # pour utiliser de la transparence dans la couleur de remplissage des points
  ) +
  facet_grid(region ~ mag_catego) +
  coord_quickmap() # puisque les points représentent des coordonnées géographiques
  # ou coord_map() si le package mapproj est installé

Question 7 :

Mapping et Facets avec ggplot2

Reproduisez avec ggplot2 l’équivalent de l’« exemple complet de représentation simultané de quatre variables » produit avec le package lattice à la section 3 des notes de cours intitulées « Graphiques en R » (https://stt4230.rbind.io/communication_resultats/graphiques_r/#package-lattice), soit le graphique suivant.

Solution :

ggplot(data = quakes) +
  geom_point(mapping = aes(x = depth, y = stations, colour = mag_catego)) +
  facet_wrap(facets = ~ region) +
  labs(
    title = "Nombres de stations ayant rapporté les séismes\nselon différentes variables",
    x = "profondeur",
    y = "nombre de stations",
    colour = "classe de\nmagnitude"
  )

Pour obtenir une mise en forme un peu plus similaire à celle de lattice :

library(ggplot2)
ggplot(data = quakes) +
  geom_point(
    mapping = aes(x = depth, y = stations, colour = mag_catego), 
    shape = 1                                # modification du type des points
  ) +
  facet_wrap(facets = ~ region) +
  labs(
    title = "Nombres de stations ayant rapporté les séismes\nselon différentes variables",
    x = "profondeur",
    y = "nombre de stations",
    colour = "classe de\nmagnitude"
  ) +
  theme_bw() +                               # thème en noir et blanc
  theme(
    plot.title = element_text(hjust = 0.5),  # permet de centrer le titre
    legend.title.align = 0.5                 # permet de centrer le nom de la légende
  ) +
  scale_colour_manual(                       # modification des couleurs des points
    values = c("4" = "dodgerblue", "5" = "magenta", "6" = "darkgreen")
  )

Question 8 :

Boxplots et mise en forme avec ggplot2

Reproduisez avec ggplot2 l’équivalent du premier exemple de graphique produit avec le système graphique R de base dans les notes de cours, soit le graphique suivant.

Solution :

ggplot(
    data = quakes,
    mapping = aes(x = mag_catego, y = stations)
  ) +
  geom_boxplot() +
  geom_jitter(
    mapping = aes(colour = region),
    alpha = 0.2, 
    width = 0.3
  ) +
  labs(
    title = "Nombres de stations ayant rapporté les séismes\nselon leurs classes de magnitude",
    x = "classe de magnitude",
    y = "nombre de stations ",
    colour = "région"
  )

Pour obtenir une mise en forme un peu plus similaire à celle du système graphique R de base :

ggplot(
    data = quakes,
    mapping = aes(x = mag_catego, y = stations)
  ) +
  geom_boxplot(
    size = 0.8,
    outlier.shape = 1,
    fill = grey(0.8)
  ) +
  geom_jitter(
    mapping = aes(colour = region),
    alpha = 0.2, 
    width = 0.3
  ) +
  labs(
    title = "Nombres de stations ayant rapporté les séismes\nselon leurs classes de magnitude",
    x = "classe de magnitude",
    y = "nombre de stations ",
    colour = "région"
  ) +
  scale_colour_manual(
    values = c("Ouest" = "red", "Est" = "blue")
  ) +
  scale_y_continuous(breaks = seq(20, 120, 20)) +  # modification des marques de graduation sur l'axe des y
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5),
    legend.position = c(0.9, 0.2),
    legend.title.align = 0.5,
    axis.text = element_text(colour = "black"),
    axis.text.y = element_text(angle = 90, hjust = 0.5)
  )

Remarque : Les limites des boîtes ne sont pas calculées exactement de la même façon par ggplot2::geom_boxplot et par graphics::boxplot. C’est mentionné ici https://ggplot2.tidyverse.org/reference/geom_boxplot.html à la section « Summary statistics ». En conséquence, les moustaches n’ont pas tout à fait les mêmes longueurs et un des deux graphiques peut contenir plus d’observations extrêmes. Ajuster le style des boxplot de ggplot2 n’est pas simple (voir https://stackoverflow.com/questions/53170465/how-to-make-a-base-r-style-boxplot-using-ggplot2 pour des pistes de solutions).

Question 9 :

Ajout d’éléments à un graphique produit avec ggplot2

Cette question est similaire à la question 4, mais nous allons utiliser des fonctions du package ggplot2 au lieu de fonctions du système graphique R de base. Aussi, le graphique sera un peu modifié pour intégrer la variable catégorique region.

Produisez avec des fonctions du package ggplot2 un diagramme de dispersion des observations de la variable stations en fonction de la variable mag, du jeu de données quakes. Les points doivent avoir une couleur différente selon la variable region.

Ensuite, ajoutez au diagramme :

  1. la courbe de régression de la variable stations en fonction de la variable mag (une seule courbe, donc ne pas tenir compte de la variable région pour cette couche);

  2. deux courbes de lissage par la méthode "loess", soit une pour chaque niveau de la variable region, accompagnées de région de confiance (avec le niveau de confiance par défaut).

  3. une légende (elle sera ajoutée automatiquement par ggplot2 si vous faites bien les choses) spécifiant la signification des différentes couleurs de points (et peut-être de lignes) et permettant de distinguer la courbe de régression des courbes de lissage;

  4. l’erreur standard résiduelle (residual standard error) de la régression, en bas à droite, avec le libellé \(\sigma_e = 11.5\) (truc : ne tentez pas d’utiliser directement l’objet outregRSE cette fois et allez voir les derniers exemples de la fiche d’aide suivante : https://ggplot2.tidyverse.org/reference/annotate.html) ?

Solution :

ggplot(
    data = quakes, 
    mapping = aes(x = mag, y = stations)
  ) +
  geom_point(
    mapping = aes(colour = region)
  ) +
  geom_smooth(  # i. Ajout de la courbe de régression
    mapping = aes(linetype = "droite de régression"),  # linetype spécifié pour obtenir légende
    method = "lm",
    se = FALSE
  ) +
  geom_smooth(  # ii. Ajout des deux courbes de lissage
    mapping = aes(colour = region, linetype = "courbe de lissage"),  # linetype spécifié pour obtenir légende
    method = "loess"
  ) +
  labs(
    x = "magnitude du séisme",
    y = "nombre de stations rapportant le séisme",
    colour = "région",
    linetype = ""
  ) +
  theme(
    legend.key.width = unit(1.5, "cm"),  # pour élargir les extraits de droites dans la légende 
  ) +
  annotate(
    geom = "text", 
    x = 6.2, y = 10, 
    label = "sigma[e] == 11.5", 
    parse = TRUE
  )
# Grâce à parse = TRUE, la valeur donnée à label est interprétée comme
# une expression plotmath (annotation mathématique).
# Différence par rapport aux annotations mathématiques avec le système de base :
# il ne faut pas envelopper la valeur d'un appel à la fonction expression,
# mais il faut l'encadrer de guillemets.


Questions plus poussées avec le système graphique de base

À partir d’ici, les questions sont plus difficiles. Faites-les seulement si votre niveau en R est déjà assez bon.

Données sur la production de fruits dans le monde

Reprenons encore une fois les données sur la production de fruits dans le monde enregistrées dans le fichier fruitsFinal.rds. Ce fichier est téléchargeable à partir du site web du cours dans le contenu de cette semaine.

Les variables de ces données sont ainsi définies :

  • Region : Région du monde, presque toujours un pays, mais parfois un territoire rattaché à un pays;
  • Item : Nom du ou des fruits, ou du groupe de fruits;
  • Annee : Année (2014, 2015 ou 2016);
  • Production : Quantité de fruits produits, en tonnes;
  • Superficie : Superficie de la région, en km2.

Téléchargez le fichier fruitsFinal.rds et chargez en R l’objet qu’il contient. Nommez cet objet fruits.

fruits <- readRDS("fruitsFinal.rds")
str(fruits)
## 'data.frame':    7758 obs. of  5 variables:
##  $ Region    : chr  "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ...
##  $ Item      : chr  "Fruits, agrumes nda" "Fruits, agrumes nda" "Figues" "Figues" ...
##  $ Annee     : int  2016 2015 2014 2015 2015 2014 2015 2014 2016 2014 ...
##  $ Production: num  17472 3043 6100 4996 50657 3000 87686 46186 17894 23561 ...
##  $ Superficie: num  652230 652230 652230 652230 652230 ...

Question 10 :

Diagramme en bâtons

Créez un diagramme à bâtons horizontaux des productions totales des items en 2016. Votre graphique final devrait être similaire à celui-ci.

Pour y arriver, suivez les étapes suivantes.

  1. Créez le jeu de données à représenter. Il faut :

    • partir du jeu de données fruits,
    • se limiter aux données pour l’année 2016,
    • calculer les productions totales pour tous les items (somme sur les valeurs pour les différentes régions),
    • ordonner le data frame obtenu en ordre décroissant de production totale.
  2. Créez un simple diagramme en bâtons horizontaux des productions totales obtenues en a.

  3. Dans le graphique, affichez les noms d’items à l’horizontale

  4. Dans le graphique, ajustez la marge de gauche pour voir les noms d’items en entier.

  5. Modifiez les données présentées dans le dernier graphique afin de regrouper tous les items avec une production totale inférieure à 1 million de tonnes en une seule catégorie. Ensuite, générez de nouveau votre graphique. Assurez-vous que le bâton pour la nouvelle catégorie se retrouve dans le bas du graphique. Ainsi, les items devraient être présentés dans le graphique dans l’ordre du plus produit (bâton du haut) au moins produit (avant dernier bâton dans le bas), suivis de la nouvelle catégorie.

  6. Dans le graphique, ajoutez un titre centré par rapport à la fenêtre graphique et non par rapport à la zone graphique.

  7. Dans le graphique, omettez tout affichage de l’axe des abscisses (l’axe des x).

  8. Dans le graphique, ajoutez à droite des bâtons les valeurs de production totale en millions de tonnes (modification la plus compliquée à réaliser).

Solution :

# Extraction de l'an 2016
fruits2016 <- subset(fruits, subset = Annee == 2016)
# Calcul des productions totales par item
fruits2016Item <- aggregate(Production ~ Item, data = fruits2016, FUN = sum, na.rm = TRUE)
# Pour ordonner les données
fruits2016Item <- fruits2016Item[order(fruits2016Item$Production),]
# Diagramme en bâtons horizontaux des productions totales pour tous les items
barplot(height = fruits2016Item$Production, 
        names.arg = fruits2016Item$Item, 
        horiz = TRUE)

# Comment mettre les noms d'items horizontaux ?
barplot(height = fruits2016Item$Production, 
        names.arg = fruits2016Item$Item, 
        horiz = TRUE, 
        las = 1)  ## grâce à l'argument las

# Comment ajuster la marge pour voir les noms d'items entiers ?
par.default <- par(mar = c(3, 20, 0, 2))  ## grâce au paramètre graphique mar
barplot(height = fruits2016Item$Production, 
        names.arg = fruits2016Item$Item, 
        horiz = TRUE, 
        las = 1)
par(par.default)

# Comment regrouper les items les moins produits en une seule catégorie ?
# Il faut créer un autre jeu de données.
fruits2016Item$ItemAutre <- ifelse(fruits2016Item$Production < 1e6, "Autres", fruits2016Item$Item)
# Calculer la production totale pour la catégorie Autres :
prodAutre <- sum(fruits2016Item$Production[fruits2016Item$ItemAutre ==  "Autres"])
# Combiner les bonnes observations :
fruits2016ItemAutre <- subset(fruits2016Item, 
                              subset = ItemAutre != "Autres", 
                              select = -Item)
fruits2016ItemAutre <- rbind(data.frame(Production = prodAutre, ItemAutre = "Autres"), fruits2016ItemAutre)
# Je n'ai pas utilisé aggregate afin de contrôler l'ordre des observations.

par.default <- par(mar = c(3,20,0,2))
barplot(height = fruits2016ItemAutre$Production, 
        names.arg = fruits2016ItemAutre$ItemAutre, 
        horiz = TRUE, 
        las = 1)
par(par.default)

# Comment ajouter un titre centré par rapport à la fenêtre
# et non par rapport à la zone graphique ?
par.default <- par(mar = c(3,20,0,2),  oma = c(0,0,2,0))
barplot(height = fruits2016ItemAutre$Production, 
        names.arg = fruits2016ItemAutre$ItemAutre, 
        horiz = TRUE, 
        las = 1)
mtext("Production mondiale de fruits en 2016 (millions de tonnes)",  ## avec mtext
      outer = TRUE, cex = 1.5)
par(par.default)

# Comment ne pas afficher l'axe des x ?
par.default <- par(mar = c(1,20,0,2),  oma = c(0,0,2,0))
barplot(height = fruits2016ItemAutre$Production, 
        names.arg = fruits2016ItemAutre$ItemAutre, 
        horiz = TRUE, 
        las = 1,
        xaxt = "n")  ## grâce à l'argument xaxt
mtext("Production mondiale de fruits en 2016 (millions de tonnes)", 
      outer = TRUE, cex = 1.5)
par(par.default)

# Comment ajouter les quantités à droite des bâtons ?
par.default <- par(mar = c(1,20,0,2), oma = c(0,0,2,0))
br <- barplot(height = fruits2016ItemAutre$Production, 
              names.arg = fruits2016ItemAutre$ItemAutre, 
              horiz = TRUE, 
              las = 1, 
              xaxt = "n", 
              xlim = c(0,max(fruits2016ItemAutre$Production)*1.3))
text(x = fruits2016ItemAutre$Production, y = br, pos = 4,        ## avec text
     labels = round(fruits2016ItemAutre$Production/1000000, 2))
mtext("Production mondiale de fruits en 2016 (millions de tonnes)", 
      outer = TRUE, cex = 1.5)
par(par.default)

Question 11 :

Représentation d’une fonction mathématique

Reproduisez le graphique suivant, qui s’est déjà trouvé sur la page d’accueil du r-project (et qui est encore disponible en ligne à l’adresse suivante : https://www.r-project.org/screenshots/power.png).

Solution :

# Création de la zone graphique et du nom des axes
plot(
  x = c(-3, 8), y = c(0, 0.4), 
  type = "n", 
  ylab = "Density",
  xlab = expression(Z == frac(mu[1] - mu[2], sigma/sqrt(n)))
)

# Zones ombragées
x <- seq(-4, 9, 0.01)
y1 <- dnorm(x)
y2 <- dnorm(x, 3)
polygon(
  x = c(x[x <= -2], -2, -3),
  y = c(y1[x <= -2], 0, 0), 
  col = "grey30", border = NA
)
polygon(
  x = c(x[x >= 2], 6, 2), 
  y = c(y2[x >= 2], 0, 0), 
  col = "grey", border = NA
)
polygon(
  x = c(x[x >= 2], 3, 2), 
  y = c(y1[x >= 2], 0, 0), 
  col = "grey30", border=NA
)

# Courbes de densités (tracées après pour voir le contour noir)
curve(expr = dnorm(x), from = -3, to = 8, n = 1001, add = TRUE)
curve(expr = dnorm(x, 3), from = -3, to = 8, n = 1001, add = TRUE)

# Texte ajouté dans le graphique
text(
  x = 0, y = 0.2, 
  labels = expression(H[0] : mu[1] == mu[2])
)
text(
  x = 3, y = 0.2, 
  labels = expression(H[1] : mu[1] == mu[2] + delta)
)

# Légende
legend(
  x = 4, y = 0.4, 
  legend = c(
    expression(P( abs(Z) > 1.96 , H[1]) == 0.85),
    expression(P( abs(Z) > 1.96 , H[0]) == 0.05)
  ),
  fill = c("grey30", "grey"), bty = "n"
)