Soit la matrice suivante :
mat <- matrix(c(2, 5, 3, 4, 6, 5, 4, 3, 1, 2, 9, 8), nrow = 3, ncol = 4)
mat## [,1] [,2] [,3] [,4]
## [1,] 2 4 4 2
## [2,] 5 6 3 9
## [3,] 3 5 1 8Pour chaque colonne de cette matrice (sauf la dernière), calculez sa somme terme à terme avec la colonne suivante (colonne 1 + colonne 2, colonne 2 + colonne 3, colonne 3 + colonne 4), de façon à obtenir le résultat suivant :
## [,1] [,2] [,3]
## [1,] 6 8 6
## [2,] 11 9 12
## [3,] 8 6 9Votre solution doit tenir en une seule commande et ne pas contenir de boucle.
Remarque : Ce calcul ressemble au calcul effectué par la fonction
diff lorsqu’elle reçoit en entrée une matrice. Cependant,
des additions sont demandées au lieu de soustractions, et les opérations
doivent être effectuées avec des paires de colonnes au lieu de paires de
lignes.
Divisez chaque élément de la matrice mat par le numéro
de sa ligne et arrondissez le résultat obtenu à l’entier supérieur, de
façon à obtenir le résultat suivant :
## [,1] [,2] [,3] [,4]
## [1,] 2 4 4 2
## [2,] 3 3 2 5
## [3,] 1 2 1 3Votre solution doit tenir en une seule commande et ne pas contenir de boucle.
ceiling(mat/(1:nrow(mat)))
## [,1] [,2] [,3] [,4]
## [1,] 2 4 4 2
## [2,] 3 3 2 5
## [3,] 1 2 1 3
# La règle de recyclage s'applique en transformant le vecteur 1:nrow(mat) en une
# matrice de la même dimension que mat. La matrice est remplie colonne par colonne
# en répétant le vecteur 1:nrow(mat) autant de fois que nécessaire. Ensuite une
# division terme à terme est effectuée entre les deux matrice. Finalement, la fonction
# ceiling effectue un arrondissement pour chacun des termes de la matrice obtenue.airqualityLes questions suivantes vous feront manipuler les données stockées
dans le data frame airquality, inclus dans le package
datasets de l’installation de base de R. Si ce package est
chargé dans votre session R (ce qui est fait par défaut lors de
l’ouverture d’une nouvelle session R), vous avez accès à l’objet
airquality. Il contient des observations de quantité
d’ozone et de radiation solaire dans l’air de New York pour chaque jour
des mois de mai à septembre en 1973. Le data frame comprend aussi des
mesures de vitesse du vent et de température.
Est-ce que le jeu de données airquality contient des
données manquantes ? Si oui, pour quelles variables ?
Bien qu’il soit possible de répondre à ces questions uniquement en regardant les données, vous devez écrire du code R qui répond à ces questions.
str(airquality)
## 'data.frame': 153 obs. of 6 variables:
## $ Ozone : int 41 36 12 18 NA 28 23 19 8 NA ...
## $ Solar.R: int 190 118 149 313 NA NA 299 99 19 194 ...
## $ Wind : num 7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
## $ Temp : int 67 72 74 62 56 66 65 59 61 69 ...
## $ Month : int 5 5 5 5 5 5 5 5 5 5 ...
## $ Day : int 1 2 3 4 5 6 7 8 9 10 ...
# Est-ce que le jeu de données airquality contient des données manquantes?
any(is.na(airquality))
## [1] TRUE
# Réponse : oui, au moins une donnée manquante est présente dans ce data frame.
# Si oui, pour quelles variables?
names(which(apply(is.na(airquality), MARGIN = 2, FUN = any)))
## [1] "Ozone" "Solar.R"
# Réponse : les variables Ozone et Solar.RSi nous voulions traiter ces données, il faudrait décider quoi faire avec les données manquantes. La solution la plus simple serait de les omettre.
Lesquelles des commandes suivantes permettent de retirer complètement
les lignes contenant au moins une donnée manquante du data frame
airquality ?
# a)
airquality[airquality$Ozone != NA & airquality$Solar.R != NA, ]
# b)
airquality[! is.na(airquality$Ozone) && ! is.na(airquality$Solar.R), ]
# c)
airquality[! is.na(airquality$Ozone) & ! is.na(airquality$Solar.R), ]
# d)
airquality[! is.na(airquality$Ozone) | ! is.na(airquality$Solar.R), ]
# e)
airquality[! (is.na(airquality$Ozone) | is.na(airquality$Solar.R)), ]
# f)
airquality[- which(is.na(airquality$Ozone) | is.na(airquality$Solar.R)), ]
# g)
subset(airquality, subset = ! is.na(Ozone) & ! is.na(Solar.R))
# h)
na.omit(airquality)# a)
test1 <- airquality[airquality$Ozone != NA & airquality$Solar.R != NA, ]
str(test1) # ne fonctionne pas,
# c'est la fonction is.na qui sert à tester si une donnée est NA
# b)
test2 <- airquality[!is.na(airquality$Ozone) && !is.na(airquality$Solar.R), ]
str(test2) # ne fonctionne pas,
# l'opérateur logique && retourne un seul logique, l'unique TRUE retourné
# par !is.na(airquality$Ozone) && !is.na(airquality$Solar.R) est recyclé
# c)
test3 <- airquality[!is.na(airquality$Ozone) & !is.na(airquality$Solar.R), ]
str(test3) # fonctionne, la condition logique est bien écrite
# d)
test4 <- airquality[!is.na(airquality$Ozone) | !is.na(airquality$Solar.R), ]
str(test4) # ne fonctionne pas,
# retire seulement les lignes pour lesquelles les deux données sont manquantes
# e)
test5 <- airquality[!(is.na(airquality$Ozone) | is.na(airquality$Solar.R)), ]
str(test5) # fonctionne, la condition logique est bien écrite
# f)
test6 <- airquality[-which(is.na(airquality$Ozone) | is.na(airquality$Solar.R)), ]
str(test6) # fonctionne, mais l'utilisation de which n'est pas nécessaire
# -which(**) a le même effet que !(**), qui est plus succinct
# g)
test7 <- subset(airquality, subset = !is.na(Ozone) & !is.na(Solar.R))
str(test7) # fonctionne, la condition logique est bien écrite
# la fonction subset permet des écritures plus courtes, car les noms des
# colonnes du data frame fourni en premier argument sont directement
# accessibles (donc plus besoin des airquality$ dans la condition logique)
# h)
test8 <- na.omit(airquality)
str(test8) # fonctionne -> solution la plus générale, car elle ne dépend pas des noms des
# variables sur lesquelles tester la présence de NA, ni du nombre de variablesSupposez que vous devez programmer vous-même le calcul de l’écart-type. La formule mathématique à programmer est celle qui se trouve dans la section « Background » de la page web suivante : https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation#Background
Contraintes :
mean ou sd par exemple), mais vous pouvez
utiliser des fonctions pour des calculs mathématiques (sum
par exemple).À faire :
Écrivez un programme R qui calcule l’écart-type des observations de
la variable Ozone dans le data frame
airquality. Pour simplifier la programmation, vous pouvez
travailler sur le vecteur x suivant :
x <- airquality$Ozonex <- airquality$Ozone
# Il faut d'abord calculer la moyenne
n <- sum(!is.na(x))
# Dénominateur = nombre d'observations non manquantes
xMoy <- sum(x, na.rm = TRUE)/n
xMoy
## [1] 42.12931
# Ensuite nous calculons l'écart-type
xSD <- sqrt(sum((x - xMoy)^2, na.rm = TRUE)/(n - 1))
xSD
## [1] 32.98788
# Est-ce que nous obtenons le même résultat que la fonction sd?
sd(x, na.rm = TRUE)
## [1] 32.98788
# oui!
# Autre solution
y <- na.omit(x)
n <- length(y)
yMoy <- sum(y)/n
ySD <- sqrt(sum((y - yMoy)^2)/(n - 1))
ySD
## [1] 32.98788Soit X une variable aléatoire suivant une distribution Poisson de
paramètre \(\lambda\) (argument
lambda) égal à 5. Combien vaut la probabilité que cette
variable aléatoire prenne la valeur \(x\), pour \(x \in
\{3, 4, 5\}\) ?
Soit X une variable aléatoire suivant une distribution exponentielle.
Combien vaut la probabilité que cette variable aléatoire prenne une
valeur inférieure ou égale à 3 si le paramètre \(\lambda\) (argument rate) de
la distribution vaut 0.2 et s’il vaut 0.5 ?
Supposons que vous venez d’effectuer un test d’hypothèses. Sous \(H_0\), l’hypothèse nulle du test, la statistique du test, notée \(X\), suit une loi de Fisher à 2 et 10 degrés de liberté. La valeur observée de la statistique de test est \(x_{obs} = 4.3\). Le seuil observée (p-value) du test est donc la probabilité suivante : \(P(X \ge 4.3)\) sachant que \(X \sim F_{2,10}\). Obtenez la valeur de ce seuil observé.
Note : Étant donné que \(X\) est une variable aléatoire continue, \(P(X = x)\) vaut approximativement 0.
Si vous choisissiez d’utiliser une valeur critique plutôt qu’un seuil observée pour émettre une conclusion dans le test d’hypothèses précédent, vous compareriez alors la valeur observée de votre statistique de test à un quantile de la loi Fisher à 2 et 10 degrés de liberté. Pour un seuil de 5%, cette valeur critique serait donc la valeur \(q\) pour laquelle \(P(X \ge q) = 0.05\). Calculez cette valeur critique.
Générez n = 100 observations tirées d’une loi normale standard.
Tracez ensuite l’histogramme de ces observations, auquel vous devez
superposer la courbe de densité de la loi normale standard. Voici le
code pour produire ce graphique, en supposant que les observations
générées sont stockées dans un vecteur nommé obs
hist(obs, freq = FALSE)
curve(dnorm, add = TRUE)Produisez ensuite le même graphique, mais à partir de n = 1000 observations tirées d’une loi normale standard.
Puis refaites encore le graphique, cette fois à partir de n = 10000 observations, toujours tirées d’une loi normale standard.
Qu’observez-vous ?
n <- 100 # nombres d'observations à générer
n <- 1000
n <- 10000
obs <- rnorm(n) # génération aléatoire
# Comparons les observations aléatoires obtenues à la distribution
# théorique en superposant une courbe de densité de la loi normale
# à l'histogramme des observations.
hist(obs, freq = FALSE)
curve(dnorm, add = TRUE)
# Nous pouvons observer l'effet de la loi des grands nombres :
# plus n est grand, plus la distribution empirique s'approche
# de la distribution théorique.Avant de générer pseudo-aléatoirement 6 observations d’une loi de Student à 15 degrés de liberté, je vais fixer le germe du générateur de nombres pseudo-aléatoires avec la commande :
set.seed(35)Quel sera le germe utilisé pour cette génération pseudo-aléatoire ?
Pour les questions suivantes, vous devez importer en R les données du fichier student-mat.csv téléchargeable sur le site web du cours. Il s’agit de résultats d’étudiants portugais dans un cours de mathématiques.
Source des données et description des variables : https://archive.ics.uci.edu/ml/datasets/Student+Performance#
student <- read.csv2("student-mat.csv")
str(student)
## 'data.frame': 395 obs. of 33 variables:
## $ school : Factor w/ 2 levels "GP","MS": 1 1 1 1 1 1 1 1 1 1 ...
## $ sex : Factor w/ 2 levels "F","M": 1 1 1 1 1 2 2 1 2 2 ...
## $ age : int 18 17 15 15 16 16 16 17 15 15 ...
## $ address : Factor w/ 2 levels "R","U": 2 2 2 2 2 2 2 2 2 2 ...
## $ famsize : Factor w/ 2 levels "GT3","LE3": 1 1 2 1 1 2 2 1 2 1 ...
## $ Pstatus : Factor w/ 2 levels "A","T": 1 2 2 2 2 2 2 1 1 2 ...
## $ Medu : int 4 1 1 4 3 4 2 4 3 3 ...
## $ Fedu : int 4 1 1 2 3 3 2 4 2 4 ...
## $ Mjob : Factor w/ 5 levels "at_home","health",..: 1 1 1 2 3 4 3 3 4 3 ...
## $ Fjob : Factor w/ 5 levels "at_home","health",..: 5 3 3 4 3 3 3 5 3 3 ...
## $ reason : Factor w/ 4 levels "course","home",..: 1 1 3 2 2 4 2 2 2 2 ...
## $ guardian : Factor w/ 3 levels "father","mother",..: 2 1 2 2 1 2 2 2 2 2 ...
## $ traveltime: int 2 1 1 1 1 1 1 2 1 1 ...
## $ studytime : int 2 2 2 3 2 2 2 2 2 2 ...
## $ failures : int 0 0 3 0 0 0 0 0 0 0 ...
## $ schoolsup : Factor w/ 2 levels "no","yes": 2 1 2 1 1 1 1 2 1 1 ...
## $ famsup : Factor w/ 2 levels "no","yes": 1 2 1 2 2 2 1 2 2 2 ...
## $ paid : Factor w/ 2 levels "no","yes": 1 1 2 2 2 2 1 1 2 2 ...
## $ activities: Factor w/ 2 levels "no","yes": 1 1 1 2 1 2 1 1 1 2 ...
## $ nursery : Factor w/ 2 levels "no","yes": 2 1 2 2 2 2 2 2 2 2 ...
## $ higher : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...
## $ internet : Factor w/ 2 levels "no","yes": 1 2 2 2 1 2 2 1 2 2 ...
## $ romantic : Factor w/ 2 levels "no","yes": 1 1 1 2 1 1 1 1 1 1 ...
## $ famrel : int 4 5 4 3 4 5 4 4 4 5 ...
## $ freetime : int 3 3 3 2 3 4 4 1 2 5 ...
## $ goout : int 4 3 2 2 2 2 4 4 2 1 ...
## $ Dalc : int 1 1 2 1 1 1 1 1 1 1 ...
## $ Walc : int 1 1 3 1 2 2 1 1 1 1 ...
## $ health : int 3 3 3 5 5 5 3 1 1 5 ...
## $ absences : int 6 4 10 2 4 10 0 6 0 0 ...
## $ G1 : int 5 5 7 15 6 15 12 6 16 14 ...
## $ G2 : int 6 5 8 14 10 15 12 5 18 15 ...
## $ G3 : int 6 6 10 15 10 15 11 6 19 15 ...Dans ces données, est-ce qu’il y a une corrélation de Pearson
significativement différente de 0 entre les résultats à la première
période (variable G1) et ceux à la deuxième période
(variable G2) ?
# 1- Visualisation des données
pairs(formula = ~ G1 + G2 + G3, data = student)
# 2- Calcul de corrélations
cor(student[, c("G1", "G2", "G3")])
# 3- Test demandé
cor.test(formula = ~ G1 + G2, data = student)
## Pearson's product-moment correlation
##
## data: G1 and G2
## t = 32.278, df = 393, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8226117 0.8770475
## sample estimates:
## cor
## 0.8521181
# Avec un p-value de 2.2e-16, oui, cette corrélation est
# significativement différente de 0.Ajustez un modèle linéaire classique (avec erreurs normales) sur ces
données. Dans ce model, la note finale (variable G3) des
étudiants est la variable réponse. Ce modèle doit comporter un terme
pour l’école fréquentée par les étudiants (variable catégorique
school) et un autre terme pour l’âge des étudiants
(variable numérique age).
# ---- 1- Visualisation des données : relations bivariées ----
# G3 versus school
boxplot(G3 ~ school, data = student)
# G3 versus age
plot(formula = G3 ~ age, data = student)
# ou
plot(x = student$age, y = student$G3)
# ou
with(student, plot(age, G3))
# Ajout d'un peu de bruit pour mieux voir les observations supperposées
with(student, plot(jitter(age), jitter(G3)))
# Autre option
boxplot(formula = G3 ~ age, data = student)
# Attention, le boxplot pour age = 20 représente en fait seulement 3 observations
# ---- 2- Ajustement du modèle ----
mod <- lm(formula = G3 ~ school + age, data = student)
mod
## Call:
## lm(formula = G3 ~ school + age, data = student)
##
## Coefficients:
## (Intercept) schoolMS age
## 20.4924 0.2661 -0.6054
# Comment visualiser les résultats des tests sur les coefficients ?
summary(mod)
## Call:
## lm(formula = G3 ~ school + age, data = student)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.4111 -1.5948 0.1943 3.1943 9.6160
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.4924 3.2018 6.400 4.43e-10 ***
## schoolMS 0.2661 0.7677 0.347 0.72904
## age -0.6054 0.1932 -3.133 0.00186 **
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 4.532 on 392 degrees of freedom
## Multiple R-squared: 0.02641, Adjusted R-squared: 0.02144
## F-statistic: 5.316 on 2 and 392 DF, p-value: 0.005273
# ---- 3- Validation des postulats du modèle ----
# Je n'entre pas dans les détails ici. Voici un exemple non complet.
plot(formula = G3 ~ age, data = student)
abline(a = coef(mod)["(Intercept)"], b = coef(mod)["age"])
plot(mod)
# ---- 4- Modification du modèle ou des données au besoin ----
# Exemples questions 11 et 12Selon le modèle ajusté à la question précédente, est-ce que la note finale des étudiants varie significativement au seuil de 5% selon l’école fréquentée ?
coef(summary(mod))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.4924417 3.2017529 6.400382 4.430048e-10
## schoolMS 0.2661383 0.7677302 0.346656 7.290357e-01
## age -0.6054218 0.1932361 -3.133067 1.859921e-03
coef(summary(mod))["schoolMS", "Pr(>|t|)"]
## [1] 0.7290357
coef(summary(mod))["schoolMS", "Pr(>|t|)"] <= 0.05
## [1] FALSE
# Seuil observé du test sur schoolMS : 0.7290357 > 0.05
# Réponse : NonSelon ce modèle, est-ce que la note finale des étudiants varie significativement au seuil de 5% selon leur âge ?
coef(summary(mod))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.4924417 3.2017529 6.400382 4.430048e-10
## schoolMS 0.2661383 0.7677302 0.346656 7.290357e-01
## age -0.6054218 0.1932361 -3.133067 1.859921e-03
coef(summary(mod))["age", "Pr(>|t|)"]
## [1] 0.001859921
coef(summary(mod))["age", "Pr(>|t|)"] <= 0.05
## [1] TRUE
# Seuil observé du test sur age : 0.00186 <= 0.05
# Réponse : OuiAjustez de nouveau le même modèle qu’à la question 8, en vous limitant cette fois aux étudiants pour lesquels :
Selon ce modèle, est-ce que la note finale des étudiants varie significativement au seuil de 5% selon leur âge ?
# - Retirer les individus ayant obtenus 0 ainsi que ceux de 20 ans ou plus
mod <- lm(
formula = G3 ~ school + age,
data = student,
subset = G3 > 0 & age < 20
)
plot(
formula = G3 ~ age,
data = student,
subset = G3 > 0 & age < 20
)
abline(
a = coef(mod)["(Intercept)"],
b = coef(mod)["age"]
)
summary(mod)
## Call:
## lm(formula = G3 ~ school + age, data = student, subset = G3 >
## 0 & age < 20)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.4143 -2.1347 -0.4143 2.2441 8.2441
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.5379 2.5465 6.887 2.66e-11 ***
## schoolMS -0.2983 0.5808 -0.513 0.6079
## age -0.3602 0.1544 -2.332 0.0202 *
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 3.184 on 349 degrees of freedom
## Multiple R-squared: 0.02129, Adjusted R-squared: 0.01568
## F-statistic: 3.795 on 2 and 349 DF, p-value: 0.02341
# Seuil observé du test sur age : 0.0202 < 0.05
# Effet de l'âge encore significatif, mais moins fort qu'auparavant.En revenant à l’utilisation de toutes les observations, ajustez de nouveau le modèle de la question 8 en modifiant une chose : la variable age doit être traitée en tant que variable catégorique (comme la variable school) plutôt que numérique.
Selon ce modèle, est-ce que la note finale des étudiants varie significativement au seuil de 5% selon leur âge ?
# - Traiter l'âge comme une variable catégorique
mod <- lm(
formula = G3 ~ school + factor(age),
data = student
)
summary(mod)
## Call:
## lm(formula = G3 ~ school + factor(age), data = student)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.2561 -1.4868 0.3099 2.9712 9.5132
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.2561 0.5007 22.481 < 2e-16 ***
## schoolMS 0.2033 0.7884 0.258 0.79664
## factor(age)16 -0.2273 0.6696 -0.339 0.73450
## factor(age)17 -1.0055 0.6854 -1.467 0.14319
## factor(age)18 -1.7693 0.7478 -2.366 0.01847 *
## factor(age)19 -3.0986 1.0706 -2.894 0.00401 **
## factor(age)20 2.6084 2.7165 0.960 0.33755
## factor(age)21 -4.4594 4.6292 -0.963 0.33599
## factor(age)22 -3.2561 4.5615 -0.714 0.47577
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## Residual standard error: 4.534 on 386 degrees of freedom
## Multiple R-squared: 0.0405, Adjusted R-squared: 0.02061
## F-statistic: 2.037 on 8 and 386 DF, p-value: 0.04124
# Afin de faire un test global sur le facteur age
anova(mod)
## Analysis of Variance Table
##
## Response: G3
## Df Sum Sq Mean Sq F value Pr(>F)
## school 1 16.8 16.759 0.8153 0.36713
## factor(age) 7 318.2 45.451 2.2110 0.03266 *
## Residuals 386 7935.0 20.557
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ou, sans passer par la fonction lm :
mod2 <- aov(
formula = G3 ~ school + factor(age),
data = student
)
summary(mod2)
## Df Sum Sq Mean Sq F value Pr(>F)
## school 1 17 16.76 0.815 0.3671
## factor(age) 7 318 45.45 2.211 0.0327 *
## Residuals 386 7935 20.56
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Seuil observé du test sur age : 0.0327 < 0.05
# Effet de l'âge significatif, mais encore moins fort qu'auparavant.La semaine dernière, nous avons importé en R des données sur la production de fruits dans le monde provenant du site web de la FAO (Organisation des Nations unies pour l’alimentation et l’agriculture, http://www.fao.org/faostat/fr/#data/QC).
J’ai nettoyé ces données. Seules les variables Region,
Item, Annee et Production ont été
conservées. Aussi, j’y ai ajouté une variable : Superficie.
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.J’ai enregistré une copie des données en format rds dans
le fichier fruitsFinal.rds. Téléchargez le fichier
fruitsFinal.rds à partir du site web du cours dans le
contenu de cette semaine et chargez en R l’objet qu’il contient. Nommez
cet objet fruits.
fruits <- readRDS("C:/coursR/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: int 17472 3043 6100 4996 50657 3000 87686 46186 17894 23561 ...
## $ Superficie: num 652230 652230 652230 652230 652230 ...Calculez le nombre d’observations par région dans le jeu de données
fruits.
Est-ce vrai de dire que le Canada a plus d’observations que la Bolivie?
# Fréquences des régions (nombre d'observations par région dans les données)
freqRegion <- table(fruits$Region)
str(freqRegion) # retourne un objet de classe table
## 'table' int [1:210(1d)] 42 72 54 66 46 12 18 18 63 45 ...
## - attr(*, "dimnames")=List of 1
## ..$ : chr [1:210] "Afghanistan" "Afrique du Sud" "Albanie" "Algerie" ...
freqRegion # difficile à visualiser dans ce sens, sortie non montrée ici
# nous transformons la table en data frame pour faciliter sa visualisation
freqRegionDF <- as.data.frame(freqRegion)
freqRegionDF
## Var1 Freq
## 1 Afghanistan 42
## 2 Afrique du Sud 72
## 3 Albanie 54
## 4 Algerie 66
## 5 Allemagne 46
## 6 Angola 12
## 7 Antigua-et-Barbuda 18
## 8 Arabie saoudite 18
## 9 Argentine 63
## 10 Armenie 45
## ...
nrow(freqRegionDF) # 210 régions
## [1] 210
# Pour trouver la réponse à la question concernant le Canada et la Bolivie,
# il faut d'abord réaliser que la chaîne de caractères représentant la Bolivie
# n'est pas simplement "Bolivie".
# Il s'agit plutôt de "Bolivie (Etat plurinational de)".
freqRegion[["Canada"]] > freqRegion[["Bolivie (Etat plurinational de)"]]
## FALSE
# Réponse : Faux, le Canada a moins d'observations que la Bolivie
# Nous aurions aussi pu utiliser la capacité de l'opérateur [[ à faire de l'identification
# à partir de noms partiels si la valeur TRUE est fournie à son argument exact
freqRegion[["Bolivie", exact = FALSE]]
# Nous verrons la semaine prochaine des fonctions qui permettraient aussi d'identifier
# de façon automatique l'observation associée à la Bolivie, par exemple :
freqRegion[[grep(pattern = "Bolivie", x = names(freqRegion))]]Quelle(s) région(s) présente(nt) le plus petit nombre d’observations ?
# Identification de la ou des régions avec la fréquence la plus faible
min(freqRegion)
# NON, nous venons d'identifier la fréquence minimale,
# mais pas la ou les régions associées à cette fréquence.
names(which.min(freqRegion))
## [1] "Bermudes"
# Attention!
# which.max et which.min retournent seulement la position du premier max ou min rencontré.
# Ici, le minimum est présent plus d'une fois.
names(freqRegion[freqRegion == min(freqRegion)])
## [1] "Bermudes" "Chine - RAS de Hong-Kong" "Erythree"
## [4] "Iles Vierges americaines" "Islande" "Lesotho"
## [7] "Liechtenstein" "Mongolie" "Nauru"
## [10] "Saint-Pierre-et-Miquelon"Combien de valeurs manquantes ou nulles (valeur de zéro) contient la
variable Production du jeu de données fruits
?
# Nombre de valeurs manquantes :
sum(is.na(fruits$Production))
## [1] 123
# ou encore en utilisant la sortie de summary
summary(fruits$Production)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0 1620 14247 334390 113542 79043138 123
# Nombre de valeurs nulles
sum(fruits$Production == 0, na.rm = TRUE)
## [1] 50
# Nombre de valeurs rencontrant une des deux conditions
sum(is.na(fruits$Production) | fruits$Production == 0)
## [1] 173Pour quelles régions le nombre d’observations diffère-t-il entre les années ?
# Calculons d'abord les fréquences croisées entre Region et Annee,
# c'est-à-dire le nombre d'observations par région pour chacune des années.
freqRegionAnnee <- xtabs( ~ Region + Annee, data = fruits)
# Équivalent à
# freqRegionAnnee <- table(Region = fruits$Region, Annee = fruits$Annee)
freqRegionAnnee
# Nous pouvons répondre à la question en observant le tableau de résultats.
# Mais ce tableau est plutôt gros! Il n'est pas affiché ici.
# Nous pouvons aussi faire un test automatiquement, ce qui est encore mieux.
# ---- Solution 1 ----
# Calculer des variances par ligne, puis identifier lignes pour lesquelles
# la variance est non nulles
# (variance non nulle => valeurs différentes entre les années)
test <- apply(freqRegionAnnee, MARGIN = 1, FUN = var) != 0
# Identification des régions pour lesquels la variance est non nulle
rownames(freqRegionAnnee)[test]
## [1] "Allemagne" "Australie" "Autriche"
## [4] "Belarus" "Chine, Taiwan Province de" "Danemark"
## [7] "Finlande" "Guam" "Hongrie"
## [10] "Irlande" "Japon" "Koweit"
## [13] "Republique de Coree" "Royaume-Uni" "Samoa americaines"
## [16] "Slovenie" "Soudan du Sud" "Tokelaou"
## [19] "Turquie"
# ou en utilisant la fonction which
names(which(test))
## même sortie que ci-dessus
# ---- Solution 2 ----
# Calculer le nombre de valeurs uniques par lignes, puis identifier les lignes
# avec un nombre de valeurs uniques supérieur à 1.
freqUnique <- apply(freqRegionAnnee, MARGIN = 1, FUN = unique)
# le résultat est une liste, avec un élément par région
# élément = vecteur des valeurs uniques
# Calculons la longueur de ce vecteur par région
nbreFreqUnique <- sapply(freqUnique, FUN = length)
# Indentification des régions avec plus d'une valeur de fréquence
names(which(nbreFreqUnique > 1))
## [1] "Allemagne" "Australie" "Autriche"
## [4] "Belarus" "Chine, Taiwan Province de" "Danemark"
## [7] "Finlande" "Guam" "Hongrie"
## [10] "Irlande" "Japon" "Koweit"
## [13] "Republique de Coree" "Royaume-Uni" "Samoa americaines"
## [16] "Slovenie" "Soudan du Sud" "Tokelaou"
## [19] "Turquie" Est-ce vrai de dire que la Superficie est manquante pour
6 régions dans le jeu de données fruits ?
summary(fruits$Superficie)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 12 43094 236800 1011222 752618 17098242 263
# Retourne le nombre de données manquantes, mais il peut y avoir
# plus d'une donnée manquante par région
regionNA <- unique(fruits$Region[is.na(fruits$Superficie)])
regionNA
## [1] "Guadeloupe" "Guam"
## [3] "Guyane francaise" "Martinique"
## [5] "Reunion" "Territoire palestinien occupe"
# Test formel
length(regionNA) == 6
## [1] TRUE
# Réponse : L'affirmation est vraie.Est-ce que le Canada produit des fruits clairement identifiés comme étant des oranges ?
any(fruits$Region == "Canada" & fruits$Item == "Oranges")
## [1] FALSE
# Réponse : NON, car il n'y a aucune ligne dans notre jeu de données combinant
# Région = Canada et Item = Oranges
# Remarque : Si nous avions obtenu la réponse TRUE, il aurait fallu s'assurer que
# les lignes combinant Région = Canada et Item = Orange contiennent une valeur de
# production non manquante et non nulle.Quelle(s) région(s) a(ont) produit le plus grand nombre d’items distincts en 2016 ?
# D'abord, limitons-nous à l'année 2016
fruits2016 <- subset(fruits, subset = Annee == 2016)
# La commande suivante fait le même travail.
fruits2016 <- fruits[fruits$Annee == 2016, ]
# Puis, retirons du jeu de données les lignes avec une production nulle ou manquante.
fruits2016p <- subset(fruits2016, subset = !is.na(Production) & Production != 0)
# La commande suivante fait le même travail.
fruits2016p <- fruits2016[!is.na(fruits2016$Production) & fruits2016$Production != 0, ]
# ---- Solution 1 ----
# Isolons les valeurs uniques d'item par région.
itemUnique <- tapply(fruits2016p$Item, INDEX = fruits2016p$Region, FUN = unique)
# Ensuite calculons le nombre d'items uniques par région.
nbreUnique <- sapply(itemUnique, length)
# Finalement, isolons la ou les régions associées au nombre le plus grand.
sum(nbreUnique == max(nbreUnique))
## [1] 1
# Ce nombre maximum est observé pour un seul pays.
# Nous pouvons donc utiliser which.max pour l'identifier
names(which.max(nbreUnique))
## [1] "Espagne"
# Réponse : L'Espagne
# ---- Solution 2 ----
# Fréquences croisées entre les régions et les items
croisees <- table(Region = fruits2016p$Region, Item = fruits2016p$Item)
# Vérification
# Dans le jeu de données fruits, est-ce possible de retrouver une région combinée à
# un item sur plus d'une ligne?
# Si oui, des fréquences de croisées sont supérieures à 1. Est-ce le cas?
table(croisees)
## croisees
## 0 1
## 5624 2527
# Non, croisees comprend uniquement des fréquences nulles et des 1.
# Alors la somme des lignes des fréquences croisées correspond au nombre d'items
# distincts produits par chaque région en 2016
nbreItems <- rowSums(croisees)
# Identification de la région avec le nombre max d'items
names(which.max(nbreItems))
## [1] "Espagne"Quelle(s) région(s) a(ont) produit la troisième plus grande variété d’items distincts en 2016 ?
# Nous pourrions ordonner en ordre décroissant les nombres d'items
nbreItemsSorted <- sort(nbreItems, decreasing = TRUE)
head(nbreItemsSorted, n = 15)
# Le troisième plus grand nombre d'items distincts est 28.
unique(nbreItemsSorted)[3]
# Les régions avec ce nombre d'items distincts sont les suivants :
names(nbreItemsSorted)[nbreItemsSorted == unique(nbreItemsSorted)[3]]Quelle(s) région(s) a(ont) produit la plus grande quantité de fruits en moyenne par année (moyenne, entre les trois années, de la production totale, peu importe l’item) ?
# Il faut d'abord calculer la production totale par région et année
prodTot <- aggregate(Production ~ Region + Annee,
data = fruits,
FUN = sum, na.rm = TRUE)
# Remarque : cette commande est équivalente
prodTot2 <- aggregate(fruits$Production,
by = list(Region = fruits$Region, Annee = fruits$Annee),
FUN = sum, na.rm = TRUE)
# Est-ce que les résultats sont vraiment les mêmes?
all.equal(prodTot, prodTot2)
## [1] "Names: 1 string mismatch"
## [2] "Attributes: < Component “row.names”: Numeric: lengths (627, 630) differ >"
## [3] "Component 1: Lengths (627, 630) differ (string compare on first 627)"
## [4] "Component 1: 544 string mismatches"
## [5] "Component 2: Numeric: lengths (627, 630) differ"
## [6] "Component 3: Numeric: lengths (627, 630) differ"
# Non! Les résultats ne sont pas de même dimension (prodTot2 a plus de lignes).
# Pourquoi?
head(prodTot)
## Region Annee Production
## 1 Afghanistan 2014 1714868
## 2 Afrique du Sud 2014 7060920
## 3 Albanie 2014 708365
## 4 Algerie 2014 5798981
## 5 Allemagne 2014 2821172
## 6 Angola 2014 4393539
head(prodTot2)
## Group.1 Group.2 x
## 1 Afghanistan 2014 1714868
## 2 Afrique du Sud 2014 7060920
## 3 Albanie 2014 708365
## 4 Algerie 2014 5798981
## 5 Allemagne 2014 2821172
## 6 Angola 2014 4393539
# À part des différences dans un nom de colonne, les 6 premières
# lignes des deux fichiers sont identiques.
# Afin d'identifier les lignes de prodTot2 et prodTot qui diffèrent,
# nous pourrions procéder comme suit :
# 1 - Concaténer en une seule chaîne de caractères toutes les données
# sur une ligne, et ce, pour les deux data frames
prodTotConcat <- apply(prodTot, MARGIN = 1, FUN = paste, collapse = "_")
prodTot2Concat <- apply(prodTot2, MARGIN = 1, FUN = paste, collapse = "_")
# 2 - Identifier les chaînes de caractères différentes
setdiff(prodTotConcat, prodTot2Concat)
## character(0)
# Interprétation de cette sortie :
# toutes les observations dans prodTot se trouvent dans prodTot2
setdiff(prodTot2Concat, prodTotConcat)
## [1] "Iles Vierges americaines_2014_ 0"
## [2] "Iles Vierges americaines_2015_ 0"
## [3] "Iles Vierges americaines_2016_ 0"
# Interprétation de cette sortie :
# Les trois observations de plus dans prodTot2Concat contiennent des
# valeurs totales de production nulles, en fait obtenues du calcul
sum(NA, na.rm = TRUE)
## [1] 0
# Cette différence s'explique par l'argument na.action présent dans
# la version de aggregate qui prend en entrée une formule, mais absent
# dans la version de aggregate qui prend les arguments x, by et FUN.
# La version de aggregate qui prend en entrée une formule applique
# par défaut la fonction na.omit sur le jeu de données fourni à
# l'argument data (mais restreint aux variables concernées).
# Nous pouvons travailler avec prodTot, parce que nous souhaitons laisser
# tomber les observations manquantes.
# Fin de la remarque
# Ensuite calculons la production totale moyenne par année.
# Nous pourrions utiliser tapply, by ou aggregate.
prodTotMoy <- tapply(prodTot$Production, INDEX = prodTot$Region, FUN = mean)
summary(prodTotMoy)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3 41071 465151 4071896 2386799 259149141
# Quelle(s) région(s) est(sont) associée(s) au maximum?
names(prodTotMoy[prodTotMoy == max(prodTotMoy)])
## [1] "Chine, continentale"
# Réponse : La ChineÀ faire seulement si vous n’êtes pas trop découragé!
Identifiez le ou les items les plus produits (en termes de production annuelle moyenne) par chacune des régions.
Est-ce vrai de dire que l’item le plus produit au Canada est la pomme ?
# Tout d'abord, retirons du jeu de données les observations
# avec une valeur de production manquante
fruitsNoNA <- subset(fruits, subset = !is.na(Production))
# Ensuite, calculons la production moyenne par année pour toutes
# les combinaisons de région et item.
prodMoy <- tapply(X = fruitsNoNA$Production,
INDEX = list(Region = fruitsNoNA$Region, Item = fruitsNoNA$Item),
FUN = mean)
# Le résultat est une matrice avec les régions en lignes et les items en colonnes.
# Identifions la production moyenne maximale par région.
prodMoyMax <- apply(prodMoy, MARGIN = 1, FUN = max, na.rm = TRUE)
# Comparons prodMoy à prodMoyMax pour identifier les items
# associés aux productions maximales par région.
itemMaxId <- apply(prodMoy == prodMoyMax, MARGIN = 1, FUN = which)
# Ici, la règle de recyclage a été appliquée pour répéter le vecteur
# prodMoyMax autant de fois qu'il y a de colonnes dans prodMoy.
# La fonction which appliquée ligne par ligne à la matrice logique
# a permis d'obtenir les indices des items associés aux productions
# moyennes maximales par région.
# En fait, on obtient un seul item à production maximal par région.
# Les maximums étaient donc uniques par région.
# Les indices dans itemMaxId font référence à l'odre des items dans
# les colonnes de prodMoy
itemMax <- colnames(prodMoy)[itemMaxId]
names(itemMax) <- names(itemMaxId)
# Le vecteur itemMax contient, pour chaque région, l'item avec la
# production annuelle moyenne maximale.
itemMax
## Afghanistan
## "Raisins"
## Afrique du Sud
## "Raisins"
## Albanie
## "Pasteques"
## Algerie
## "Pasteques"
## Allemagne
## "Raisins"
## Angola
## "Bananes"
## ...
itemMax[["Canada"]]
## [1] "Pommes"
# Réponse : L'affirmation est vraie.Matière non évaluée, à faire seulement si c’est pertinent pour vous.
Utilisez une fonction présentée dans le section « Algèbre linéaire » des notes de cours de cette semaine pour calculer en une seule commande la distance euclidienne entre toutes les paires d’éléments du vecteur suivant :
x <- c(5, 3, 8, 4, 6, 2)Rappel : la distance euclidienne entre deux points x1 et x2 dans un espace à une seule dimension est la racine carrée de (x2 - x1)^2 (en 1D la distance euclidienne revient donc à une différence en valeur absolue).
Vous devriez arriver aux mêmes valeurs (pas au même format de sortie) que celles obtenues avec la commande suivante :
dist(x)x <- c(5, 3, 8, 4, 6, 2)
sqrt(outer(x, x, FUN = "-")^2)
#ou
abs(outer(x, x, FUN = "-"))
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 2 3 1 1 3
## [2,] 2 0 5 1 3 1
## [3,] 3 5 0 4 2 6
## [4,] 1 1 4 0 2 2
## [5,] 1 3 2 2 0 4
## [6,] 3 1 6 2 4 0
dist(x, diag = TRUE)
## 1 2 3 4 5 6
## 1 0
## 2 2 0
## 3 3 5 0
## 4 1 1 4 0
## 5 1 3 2 2 0
## 6 3 1 6 2 4 0
# La fonction dist affiche seulement un triangle
# (le triangle inférieur par défaut) de la matrice
# de distance obtenue puisque celle-ci est symétrique.La matrice suivante est orthogonale :
A <- matrix(c(cos(2), sin(2), -sin(2), cos(2)), nrow = 2, ncol = 2)
A
## [,1] [,2]
## [1,] -0.4161468 -0.9092974
## [2,] 0.9092974 -0.4161468Une matrice orthogonale \(Q\) possède la propriété que \(Q^t Q = QQ^t = I\), où \(I\) est la matrice identité et l’exposant \(t\) représente la transposée de la matrice.
En conséquence, l’inverse de la matrice est égal à sa transposée (\(Q^{-1} = Q^t\)) pour une matrice othogonale.
Vérifiez si la matrice A respecte bien ces propriétés.