quakesPlusieurs 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")
)
plotQue produit la fonction plot lorsqu’elle reçoit en
entrée comme premier argument
ts, tel que l’objet suivantWWWusage
table à une ou deux dimensions, tel
que les objets suivantsfreq1 <- table(quakes$mag_catego)
freq2 <- table(quakes$mag_catego, quakes$region)
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ïqueAvec 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.
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"
)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
ajoutez un peu de bruit aux observations (sur les deux axes), à
l’aide de la fonction jitter.
produisez le diagramme de dispersion avec la fonction
smoothScatter,
# 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"
)Comment ajouter au diagramme de dispersion initial de la question précédente :
outreg <- lm(stations ~ mag, data = quakes)
outlowess <- lowess(quakes$mag, quakes$stations)
une légende, en haut à gauche, pour identifier les courbes ajoutées;
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
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)
)
)Reproduisez le graphique suivant, composé de plusieurs
sous-graphiques, représentant des variables du jeu de données
quakes.
# 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))ggplot2Chargement du package ggplot2 :
library(ggplot2)
ggplot2Reproduisez 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.
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éggplot2Reproduisez 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.
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")
)ggplot2Reproduisez 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.
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).
ggplot2Cette 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 :
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);
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).
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;
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)
?
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.À partir d’ici, les questions sont plus difficiles. Faites-les seulement si votre niveau en R est déjà assez bon.
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 ...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.
Créez le jeu de données à représenter. Il faut :
fruits,Créez un simple diagramme en bâtons horizontaux des productions totales obtenues en a.
Dans le graphique, affichez les noms d’items à l’horizontale
Dans le graphique, ajustez la marge de gauche pour voir les noms d’items en entier.
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.
Dans le graphique, ajoutez un titre centré par rapport à la fenêtre graphique et non par rapport à la zone graphique.
Dans le graphique, omettez tout affichage de l’axe des abscisses (l’axe des x).
Dans le graphique, ajoutez à droite des bâtons les valeurs de production totale en millions de tonnes (modification la plus compliquée à réaliser).
# 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)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).
# 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"
)