Déroulement de l’examen

Vous travaillerez de 9h15 à 12h30 au plus tard. Vous aurez accès au cours, à l’internet… La seule exigence est que vous n’utilisiez pas d’aide extérieure.

Vous devrez répondre aux questions dans ce fichier en ajoutant du texte et du code. Les questions sont sous la forme de puces :

  • Ceci est une question.

Utilisez le format “Blockquote” pour le texte de vos réponses (en mode visuel, utilisez le menu “Format” ; en mode source, commencez chaque paragraphe par > et sautez une ligne avant vos réponses). Ajoutez le code dans des bouts de code standard:

# Code
library("tidyverse")
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Tricotez souvent pour ne pas accumuler des erreurs de code.

Vous devez répondre aussi clairement que possible, en rédigeant vos réponses, mais soyez factuels et concis : inutile de montrer que vous connaissez le cours puisque vous avez accès aux documents. Votre code doit être facile à lire, avec des commentaires quand c’est utile, formaté correctement (attention aux espaces, aux indentations, à la clarté des noms des objets, etc.).

A la fin de votre travail, envoyez le fichier examen.Rmd à eric.marcon@agroparistech.fr

Les données nécessaires se trouvent dans le sous-dossier data.

L’examen comporte plusieurs questions à traiter :

  • Une régression,
  • Une AFC,
  • Une bagarre avec les données,
  • Un problème de probabilités avec du code R.

Régression

Les brochets sont des prédateurs supérieurs qui cumulent des pesticides au cours de leur vie. Vous devez quantifier le lien entre la concentration en DDT et l’âge des poissons.

Les données (âge et concentration en DDT) de 15 brochets sont dans le fichier data/Brochets.csv au format anglo-saxon.

  • Visualisez les données.
library("readr") # Déjà chargé par library("tidiverse")
read_csv("data/Brochet.csv") %>% 
  print() -> 
  brochets
## Rows: 15 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): Obs, Age, TxDDT
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## # A tibble: 15 × 3
##      Obs   Age TxDDT
##    <dbl> <dbl> <dbl>
##  1     1     2  0.2 
##  2     2     2  0.25
##  3     3     2  0.18
##  4     4     3  0.19
##  5     5     3  0.29
##  6     6     3  0.28
##  7     7     4  0.31
##  8     8     4  0.33
##  9     9     4  0.36
## 10    10     5  0.71
## 11    11     5  0.38
## 12    12     5  0.47
## 13    13     6  1.1 
## 14    14     6  0.87
## 15    15     6  0.83

La colonne Obs n’a aucune utilité.

  • Estimez un modèle linéaire simple, vérifiez les hypothèses.

Représentation grahphique :

brochets %>% 
  ggplot(aes(x = Age, y = TxDDT)) +
  geom_point() +
  geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

Estimation du modèle :

brochets_lm <- lm(TxDDT ~ Age, data = brochets)
summary(brochets_lm)
## 
## Call:
## lm(formula = TxDDT ~ Age, data = brochets)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.24133 -0.10500  0.01133  0.08300  0.30733 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.23533    0.11269  -2.088    0.057 .  
## Age          0.17133    0.02656   6.450 2.16e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1455 on 13 degrees of freedom
## Multiple R-squared:  0.7619, Adjusted R-squared:  0.7436 
## F-statistic: 41.61 on 1 and 13 DF,  p-value: 2.165e-05

Le modèle explique 76% de la variabilité. La variable age est significative. Comme la constante n’est pas différente de 0, le modèle pourrait être estimé sans elle, d’autant plus que la concentration en DDT à la naissance des brochets devrait être nulle.

Vérification de la distribution des résidus :

plot(brochets_lm, which = 1)

Les résidus ont une forme de cuvette : les points centraux ont des résidus négatifs, les extrêmes des valeurs positives.

Vérification de la normalité des résidus :

plot(brochets_lm, which = 2)

Le test de Shapiro ne rejette pas la normalité des résidus. Remarque : le nombre de points est très petit, la puissance du test est donc faible.

brochets_lm %>% residuals() %>% shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.96554, p-value = 0.7874

Effet de levier : le point 13 (6 ans, taux = 1,1), tout en haut à droite du nuage de points, a un fort effet de levier positif, compensé par celui, négatif du point 11 (le plus bas des brochets de 5 ans). Ces deux points apparaissent dans toutes les figures de diagnostic.

plot(brochets_lm, which = 5)

Conclusion : la régression ne viole pas les hypothèses de façon flagrante mais souffre de plusieurs problèmes : les résidus ne sont pas répartis de façon homogène, deux points ont un effet de levier excessif. Le fait que les ages soient des valeurs discrètes limite la pertinence du modèle.

  • Estimez le modèle qui relie le logarithme de la concentration en DDT à l’âge. Comparez-le au précédent.

Représentation graphique :

brochets %>% 
  ggplot(aes(x = Age, y = TxDDT)) +
  geom_point() +
  scale_y_log10() +
  geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

Estimation du modèle :

brochets_lm_log <- lm(log(TxDDT) ~ Age, data = brochets)
summary(brochets_lm_log)
## 
## Call:
## lm(formula = log(TxDDT) ~ Age, data = brochets)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37122 -0.15103  0.04114  0.09496  0.32278 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.44086    0.17135 -14.245 2.61e-09 ***
## Age          0.36890    0.04039   9.134 5.09e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2212 on 13 degrees of freedom
## Multiple R-squared:  0.8652, Adjusted R-squared:  0.8548 
## F-statistic: 83.43 on 1 and 13 DF,  p-value: 5.092e-07

Le modèle linéaire s’ajuste mieux après la transformation de la variable expliquée avec un R² de 87%. La constante est maintenant significative à cause de la transformation logarithmique.

Vérification de la distribution des résidus : sans amélioration par rapport au modèle original.

plot(brochets_lm_log, which = 1)

Vérification de la normalité des résidus : similaire au précédent.

plot(brochets_lm_log, which = 2)

Le test de Shapiro ne rejette pas la normalité des résidus.

brochets_lm_log %>% residuals() %>% shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.95807, p-value = 0.6589

Effet de levier : les point 13 et 11 ont toujours les mêmes effets, mais apparemment moins importants que dans le modèle original.

plot(brochets_lm_log, which = 5)

Conclusion : la régression du logarithme de la concentration en DDT fonctionne mieux que le modèle original.

  • Essayez d’améliorer la régression en permettant une relation non-linéaire par un modèle polynomial de degré 2.

Estimation du modèle :

brochets_lm_log_poly <- lm(log(TxDDT) ~ poly(Age, degree = 2), data = brochets)
summary(brochets_lm_log_poly)
## 
## Call:
## lm(formula = log(TxDDT) ~ poly(Age, degree = 2), data = brochets)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30079 -0.09397 -0.04723  0.14917  0.32431 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -0.96527    0.04875 -19.800 1.57e-10 ***
## poly(Age, degree = 2)1  2.02054    0.18881  10.701 1.71e-07 ***
## poly(Age, degree = 2)2  0.45644    0.18881   2.417   0.0325 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1888 on 12 degrees of freedom
## Multiple R-squared:  0.9093, Adjusted R-squared:  0.8942 
## F-statistic: 60.18 on 2 and 12 DF,  p-value: 5.553e-07

Le coefficient de l’âge au carré est significatif (p-value = 3%) et le R² est amélioré. Le R² ajusté, qui prend en compte l’ajout d’un paramètre, est aussi amélioré. Le modèle est donc meilleur que le précédent.

Vérification de la distribution des résidus : La distribution des résidus est beaucoup plus homogène.

plot(brochets_lm_log_poly, which = 1)

Vérification de la normalité des résidus : similaire au précédent.

plot(brochets_lm_log_poly, which = 2)

Le test de Shapiro ne rejette pas la normalité des résidus.

brochets_lm_log_poly %>% residuals() %>% shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.96903, p-value = 0.8434

Effet de levier : le point 11 ne pose plus de problème.

plot(brochets_lm_log_poly, which = 5)

Conclusion : la régression du logarithme de la concentration en DDT par l’age et l’age au carré est le meilleur modèle. Le polynôme de degré 2 a réglé le problème d’hétérogénéité de la distribution des résidus (négatifs pour les valeurs intermédiaires) en permettant d’ajuster la courbe de régression.

AFC

Les abondances de 6 espèces d’oiseaux dans 4 régions ont été tirées de Blondel et Farré (1988). Elles se trouvent dans le fichier data/Oiseaux.csv. Les fonctions que vous allez utiliser ont besoin d’un tableau avec des noms de ligne. Vous devez donc préparer les données :

read_csv2("data/Oiseaux.csv") %>% 
  print() -> 
  oiseaux
## ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.
## Rows: 6 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (1): Espèce
## dbl (4): Pologne, Bourgogne, Provence, Corse
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## # A tibble: 6 × 5
##   Espèce Pologne Bourgogne Provence Corse
##   <chr>    <dbl>     <dbl>    <dbl> <dbl>
## 1 S1          33       507       77   278
## 2 S2         157       512      164   288
## 3 S3         751       444      194   406
## 4 S4         621       221      357   494
## 5 S5         776       375      347   512
## 6 S6         854       536      580   598
  • mettez de côté dans un vecteur les noms des espèces,
especes <- oiseaux[, 1]
  • éliminez la première colonne du dataframe,
oiseaux <- oiseaux[, -1]
  • nommez les lignes du dataframe sur le modèle : rownames(nom_du_tableau) <- vecteur_des_noms.
row.names(oiseaux) <- especes$Espèce
## Warning: Setting row names on a tibble is deprecated.

Vous pouvez maintenant analyser les données.

oiseaux
## # A tibble: 6 × 4
##   Pologne Bourgogne Provence Corse
## *   <dbl>     <dbl>    <dbl> <dbl>
## 1      33       507       77   278
## 2     157       512      164   288
## 3     751       444      194   406
## 4     621       221      357   494
## 5     776       375      347   512
## 6     854       536      580   598
  • Montrez par un test statistique que les espèces d’oiseaux ne sont pas distribuées indépendamment des régions.

$ Le test du \(\chi^2\) est approprié:

chisq.test(oiseaux)
## 
##  Pearson's Chi-squared test
## 
## data:  oiseaux
## X-squared = 1286.4, df = 15, p-value < 2.2e-16

La probabilité de se tromper en rejetant l’hypothèse d’indépendance des lignes et des colonnes est très proche de 0 : les oiseaux ne sont pas distribués indépendamment des régions.

Attention : les données sont présentées sous la forme d’un tableau de contingence dont la première colonne (l’espèce) doit être retirée parce qu’elle ne contient pas de données utiles.

  • Faites une figure en mosaïque qui vous permettra de comprendre la distribution. Quelles espèces sont caractéristiques de quelles régions ?
mosaicplot(oiseaux, shade = TRUE, main = "")

Des patrons de distribution sont visibles :

  • Les espèces 1 et 2 se trouvent en Bourgogne et pas en Pologne,
  • Les espèces 3 et 5 sont en Pologne,
  • Les espèces 4 et 6 sont en Provence.
  • Effectuez une analyse factorielle des correspondances et interprétez les résultats. Vérifiez leur cohérence avec ceux de la figure en mosaïque.
library("FactoMineR")
oiseaux %>% 
  CA() ->
  oiseaux_CA

Aide : vous utiliserez la fonction CA() du package FactoMineR.

Valeurs propres :

library("factoextra")
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
oiseaux_CA %>% 
  fviz_eig()
## Warning in geom_bar(stat = "identity", fill = barfill, color = barcolor, :
## Ignoring empty aesthetic: `width`.

Le premier axe représente presque toute la variabilité. On en affichera deux pour la lisibilité.

oiseaux_CA %>% 
  fviz_ca_biplot()

L’axe 1 associe la Bourgogne aux espèces 1 et 2. Les espèces 3 à 6 sont à l’opposé. L’axe 2 montre l’opposition entre la Pologne (espèces 3 et 5) et la Provence (espèces 4 et 6). La Corse se trouve au centre de gravité du graphique, sans espèce caractéristique claire. Ces résultats sont cohérents avec ceux du graphique en mosaïque, même si les espèces 1 et 4 présentes en Corse ne sont pas visibles dans l’AFC.

On pourrait s’attendre à un gradient géographique (nord-sud, plus ou moins continental, est-ouest) mais l’AFC ne montre rien de tel. La distribution des espèces d’oiseau répond à d’autres facteurs, inexplicables avec ces données.

Bagarre

Faites une carte de la parcelle 6 de Paracou (data/Paracou6.csv) qui montrera l’emplacement du plus gros arbre de chaque famille botanique.

Aide: vous devrez d’abord trouver la taille du plus gros arbre de chaque famille, mais vous perdrez toutes les informations individuelles sur les arbres. Vous devrez donc utiliser une jointure pour retrouver les coordonnées de l’arbre. Réfléchissez bien aux colonnes qui vont servir à la jointure et nommez la bonne colonne avec le bon nom au bon moment.

Lecture des données

read_csv2("data/Paracou6.csv") ->
  paracou6
## ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.
## Rows: 3541 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (3): Family, Genus, Species
## dbl (5): SubPlot, idTree, Xfield, Yfield, CircCorr
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Sélection de l’arbre le plus gros

paracou6 %>% 
  group_by(Family) %>% 
  # Circonférence max. 
  # La colonne calculée garde le nom de la variable pour la jointure à suivre
  summarise(CircCorr = max(CircCorr)) ->
  paracou6_max

Jointure avec la table des arbres et carte

paracou6 %>% 
  inner_join(paracou6_max) %>% 
  ggplot(aes(x = Xfield, y = Yfield, color = Family, size = CircCorr)) +
  geom_point() +
  coord_fixed() +
  labs(x = "", y = "", color = "Famille", size = "Circonférence") +
  # Suppression de la légende, trop grande
  theme(legend.position = "none")
## Joining with `by = join_by(Family, CircCorr)`

Le graphique contient trop d’entrées de légende pour que cette information soit lisible : on peut supprimer la légende. L’information sur les familles ne peut pas être représentée de cette façon.

Probabilités

La méthode de Monte-Carlo permet de tester des hypothèses ou d’estimer des valeurs par simulation. Vous allez l’appliquer ici.

Le graphique suivant montre un quart de cercle de rayon 1 dans le carré de côté 1, les deux tracés à partir de l’origine du repère :

# La fonction circle calcule les coordonnées d'un point du cercle de rayon 1
# theta est l'angle de la coordonnée polaire du point
circle <- function(theta) {
  z <- exp(1i * theta)
  return(c(Re(z), Im(z)))
}
# Pour tous les angles de 0 à pi/2, par valeur de 0,01
seq(from = 0, to = pi / 2, by = 1/100) %>% 
  # Calculer les coordonnées du point sur le cercle
  sapply(FUN = circle) %>% 
  # Transposer la matrice obtenue (ligne 1 = x, ligne 2 = y)
  t() %>% 
  # La transformer en tibble, nommer les colonnes
  as_tibble(.name_repair = ~c("x", "y")) %>% 
  # Tracer le cercle
  ggplot(aes(x = x, y = y)) +
    geom_line() +
    coord_fixed() +
    # Ajouter le carré
    geom_hline(yintercept = 0:1) +
    geom_vline(xintercept = 0:1) ->
  # Enregistrer le graphique
  gg_circle

# Afficher le graphique
gg_circle

Ne modifiez pas le code ci-dessus : il produit la figure de base sur laquelle vous allez ajoutez des éléments.

  • Tirez 10 points aléatoirement (uniformément) dans le carré de côté 1 et affichez-les sur le graphique.

Aide : créez un tibble appelé mes_points dont les colonnes s’appellent x et y, dont le contenu sera un vecteur de nombres aléatoires tirés dans la loi uniforme. Ajoutez une géométrie au graphique gg_circle en précisant ses données (data = mes_points) et son esthétique (, mapping = aes(...)).

# Nombre de simulations
simulations_n <- 10
# Tirage des points
mes_points <- tibble(
  x = runif(simulations_n),
  y = runif(simulations_n)
)
# Ajout au graphique
gg_circle +
  geom_point(
    data = mes_points,
    aes(x = x, y = y)
  )

  • La probabilité qu’un point se trouve dans le quart de cercle centré sur l’origine du repère est \(\pi / 4\). Expliquez pourquoi par la géométrie.

Les points sont tirés dans le carré de surface 1. La surface du quart de cercle de rayon 1 est \(\pi / 4\). La probabilité est le rapport des surfaces.

Aide : le rayon du cercle est 1.

  • On considère être un succès qu’un point tiré se trouve dans le quart de cercle. Quelle loi suit le nombre de succès pour \(n_{sim}\) tirages ? Quelle est son espérance ?

Le nombre de succès suit une loi binomiale d’espérance \(n_{sim} * \pi / 4\).

  • On note \(n_+\) le nombre de succès. Ecrire la valeur de \(\pi\) en fonction de \(n_{sim}\) et \(n_+\)

Le nombre de succès observé est l’estimateur de l’espérance : \(n_{sim} * \pi / 4\), ce qui équivaut à \(\pi = 4 * n_+ / n_{sim}\).

  • Dans le tibble mes_points, ajoutez une colonne qui indique si chaque tirage est un succès.
mes_points %>% 
  mutate(
    # Distance du point à l'origine
    distance = sqrt(x^2 + y^2),
    # Le point est dans le cercle si distance < rayon
    succes = (distance <= 1)
  ) ->
mes_points

Aide : calculez d’abord la distance du point à l’origine puis testez si elle est inférieure au rayon du cercle (il vous faudra deux colonnes).

  • Comptez le nombre de succès pour connaître \(n_+\) et estimez \(\pi\). Comparez avec sa valeur : 3.1415927.
# Estimation de pi
4 * sum(mes_points$succes) / simulations_n
## [1] 3.6

La première décimale de \(\pi\) est correcte.

  • Tirez maintenant un million de points au lieu de 10 et donnez votre nouvelle estimation de \(\pi\).

Code complet:

# Nombre de simulations
simulations_n <- 1E6
# Tirage des points
tibble(
  x = runif(simulations_n),
  y = runif(simulations_n)
) %>% 
  mutate(
    # Distance du point à l'origine
    distance = sqrt(x^2 + y^2),
    # Le point est dans le cercle si distance < rayon
    succes = (distance <= 1)
  ) ->
  mes_points
# Estimation de pi
4 * sum(mes_points$succes) / simulations_n
## [1] 3.139592

L’estimation est plus proche de la vraie valeur mais reste assez grossière. Il faut augmenter encore le nombre de tirages, mais le temps de calcul et le besoin en mémoire pour stocker mes_points va rapidement poser problème.

Références

J. Blondel & H. Farré. The convergent trajectories of bird communities along ecological successions in european forests. Öcologia (Berlin), 75 :83–93, 1988.