Science des données biologiques

Réalisé par le service d'Écologie numérique des Milieux aquatiques, Université de Mons (Belgique)

Préambule

Si vous n’avez jamais utilisé de tutoriel “learnr”, familiarisez-vous d’abord avec son interface ici.

Conformément au RGPD (Règlement Général sur la Protection des Données), nous sommes tenus de vous informer de ce que vos résultats seront collecté afin de suivre votre progression. Les données seront enregistrées au nom de l’utilisateur apparaissant en haut de cette page. Corrigez si nécessaire ! En utilisant ce tutoriel, vous marquez expressément votre accord pour que ces données puissent être collectées par vos enseignants et utilisées pour vous aider et vous évaluer. Après avoir été anonymisées, ces données pourront également servir à des études globales dans un cadre scientifique et/ou éducatif uniquement.

Objectifs

  • Pouvoir comparer plus de deux populations simultanément en utilisant des techniques de décomposition de la variance

  • Découvrir le modèle linéaire, anciennement analyse de variance (ANOVA)

  • Savoir effectuer des tests de comparaison multiples

Distribution F

  • Représentez graphiquement la distribution F dont :
    • les degrés de liberté inter groupe est de 4
    • les degrés de liberté intra-groupe est de 96
# Fisher-Snedecor's F distribution (density probability) with parameter:
.df1 <- NUMERATOR_DF; .df2 <- DENOMINATOR_DF # numerator (.df1) and denominator (.df2) df
.col <- 1; .add <- FALSE # Plot parameters
.x <- seq(0, qf(0.999, df1 = .df1, df2 = .df2), l = 1000)  # Quantiles
.d <- function (x) df(x, df1 = .df1, df2 = .df2)           # Distribution function
.q <- function (p) qf(p, df1 = .df1, df2 = .df2)           # Quantile for lower-tail prob
.label <- bquote(F(.(.df1), .(.df2)))                      # Curve parameters
curve(.d(x), xlim = range(.x), xaxs = "i", n = 1000, col = .col,
  add = .add, xlab = "Quantiles", ylab = "Probability density") # Curve
abline(h = 0, col = "gray") # Baseline
# Fisher-Snedecor's F distribution (density probability) with parameter:
.df1 <- 4; .df2 <- 96 # numerator (.df1) and denominator (.df2) df
.col <- 1; .add <- FALSE # Plot parameters
.x <- seq(0, qf(0.999, df1 = .df1, df2 = .df2), l = 1000)  # Quantiles
.d <- function (x) df(x, df1 = .df1, df2 = .df2)           # Distribution function
.q <- function (p) qf(p, df1 = .df1, df2 = .df2)           # Quantile for lower-tail prob
.label <- bquote(F(.(.df1), .(.df2)))                      # Curve parameters
curve(.d(x), xlim = range(.x), xaxs = "i", n = 1000, col = .col,
  add = .add, xlab = "Quantiles", ylab = "Probability density") # Curve
abline(h = 0, col = "gray") # Baseline
#TODO
  • Représentez graphiquement la distribution F correspondant à une étude ayant 50 observations et 5 groupes
# Fisher-Snedecor's F distribution (density probability) with parameter:
.df1 <- 4; .df2 <- 45 # numerator (.df1) and denominator (.df2) df
.col <- 1; .add <- FALSE # Plot parameters
.x <- seq(0, qf(0.999, df1 = .df1, df2 = .df2), l = 1000)  # Quantiles
.d <- function (x) df(x, df1 = .df1, df2 = .df2)           # Distribution function
.q <- function (p) qf(p, df1 = .df1, df2 = .df2)           # Quantile for lower-tail prob
.label <- bquote(F(.(.df1), .(.df2)))                      # Curve parameters
curve(.d(x), xlim = range(.x), xaxs = "i", n = 1000, col = .col,
  add = .add, xlab = "Quantiles", ylab = "Probability density") # Curve
abline(h = 0, col = "gray") # Baseline
#TODO
  • Au seuil \(\alpha\) de 5%, Faut il rejeter \(H_0\) avec :
    • F = 2.1
    • 3 groupes
    • 55 observations

Représentez la distribution de F

# représenter graphiquement la distribution F
# 
# représenter graphiquement la distribution F
# Fisher-Snedecor's F distribution (density probability) with parameter:
.df1 <- 2; .df2 <- 52 # numerator (.df1) and denominator (.df2) df
.col <- 1; .add <- FALSE # Plot parameters
.x <- seq(0, qf(0.999, df1 = .df1, df2 = .df2), l = 1000)  # Quantiles
.d <- function (x) df(x, df1 = .df1, df2 = .df2)           # Distribution function
.q <- function (p) qf(p, df1 = .df1, df2 = .df2)           # Quantile for lower-tail prob
.label <- bquote(F(.(.df1), .(.df2)))                      # Curve parameters
curve(.d(x), xlim = range(.x), xaxs = "i", n = 1000, col = .col,
  add = .add, xlab = "Quantiles", ylab = "Probability density") # Curve
abline(h = 0, col = "gray") # Baseline
#TODO

Calculez le quantile correspondant à une aire de 5% à droite via la fonction suivantes

qf(PROBABILITES, df1 = NUMERATOR_DF, df2 =  DENOMINATOR_DF, lower.tail = TRUE)
qf(0.05, df1 = 2, df2 = 52, lower.tail = FALSE)
#TODO

Ajoutez la valeur obtenue sur la graphique avec la fonction abline()

# Fisher-Snedecor's F distribution (density probability) with parameter:
.df1 <- NUMERATOR_DF; .df2 <- DENOMINATOR_DF # numerator (.df1) and denominator (.df2) df
.col <- 1; .add <- FALSE # Plot parameters
.x <- seq(0, qf(0.999, df1 = .df1, df2 = .df2), l = 1000)  # Quantiles
.d <- function (x) df(x, df1 = .df1, df2 = .df2)           # Distribution function
.q <- function (p) qf(p, df1 = .df1, df2 = .df2)           # Quantile for lower-tail prob
.label <- bquote(F(.(.df1), .(.df2)))                      # Curve parameters
curve(.d(x), xlim = range(.x), xaxs = "i", n = 1000, col = .col,
  add = .add, xlab = "Quantiles", ylab = "Probability density") # Curve
abline(h = 0, col = "gray") # Baseline
#
# ajout du quantile correspondant à une aire de 5%
abline(v= QUANTILES, col = "red")
# représenter graphiquement la distribution F
# 
# ajout du quantile correspondant à une aire de 5%
# 
# représenter graphiquement la distribution F
# Fisher-Snedecor's F distribution (density probability) with parameter:
.df1 <- 2; .df2 <- 52 # numerator (.df1) and denominator (.df2) df
.col <- 1; .add <- FALSE # Plot parameters
.x <- seq(0, qf(0.999, df1 = .df1, df2 = .df2), l = 1000)  # Quantiles
.d <- function (x) df(x, df1 = .df1, df2 = .df2)           # Distribution function
.q <- function (p) qf(p, df1 = .df1, df2 = .df2)           # Quantile for lower-tail prob
.label <- bquote(F(.(.df1), .(.df2)))                      # Curve parameters
curve(.d(x), xlim = range(.x), xaxs = "i", n = 1000, col = .col,
  add = .add, xlab = "Quantiles", ylab = "Probability density") # Curve
abline(h = 0, col = "gray") # Baseline
# ajout du quantile correspondant à une aire de 5%
abline(v= 3.175, col = "red")
#TODO
  • Au seuil \(\alpha\) de 5%, Faut il rejeter \(H_0\) avec :
    • F = 14.4
    • 4 groupes
    • 40 observations

Calculez l’aire à droite correspondant à un quantile de 14.4

pf(QUANTILES, df1 = NUMERATOR_DF, df2 = DENOMINATOR_DF, lower.tail = TRUE)
pf(14.4, df1 = 3, df2 = 36, lower.tail = FALSE)
#TODO

Cas d’étude : La croissance des dents de cochons d’Inde

Vous allez réaliser une analyse complète sur le jeu de données portant sur la croissance des dents de cochons d’Inde. La question est : y a t’il une différence de croissance des dents en fonction de la concentration en supplément administré.

# importation
(toothgrowth <- read("ToothGrowth", package = "datasets", lang = "fr"))

Ce jeu de données comprend 3 variables :

  • len : Longueur des dents
  • supp : Supplément administré : OJ ou VC
  • dose : Dose administrée : 0.5 < 1 < 2

Vérification de l’encodage des variables

La première étape est de se demander si les variables sont convenablement encodées. La fonction glimpse() permet d’obtenir des informations sur le tableau de données.

glimpse(toothgrowth)
Observations: 60
Variables: 3
$ len  <dbl> 4.2, 11.5, 7.3, 5.8, 6.4, 10.0, 11.2, 11.2, 5.2, 7.0, 16.5,…
$ supp <fct> VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC, VC,…
$ dose <dbl> 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 1.0,…

Vous pouvez observer que la variable dose est encodée en variable numérique et non en variable facteur.

  • Réencodez la variable dose afin d’avoir une variable facteur ordonnée.

  • Visualisez à nouveau vos données avec la fonction glimpse().

[1] "as.ordered()" "dose"         "glimpse()"    "toothgrowth" 
[5] "$"           
toothgrowth <- read("ToothGrowth", package = "datasets", lang = "fr")
# réordonner la variable dose
# 
# visualiser le jeu de données
# 
# réordonner la variable dose
toothgrowth$dose <- as.ordered(toothgrowth$dose)
# visualiser le jeu de données
glimpse(toothgrowth)
#TODO

Vous pouvez également ajouter les labels et les unités aux variables de votre jeu de données avec la fonction labelise().

Variable label unité
len Longueur des dents mm
supp Supplémentation NA
dose Dose mg/J
toothgrowth <- read("ToothGrowth", package = "datasets", lang = "fr")
toothgrowth$dose <- as.ordered(toothgrowth$dose)
# Compléter la fonction
toothgrowth <- labelise(toothgrowth, self = FALSE,
  label = list(
    
  ),
  units = list(
    
  )
)
toothgrowth <- labelise(toothgrowth, self = FALSE,
  label = list(
    len = "Longueur des dents",
    supp = "Supplémentation",
    dose = "Dose"
  ),
  units = list(
    len = "mm",
    supp = NA,
    dose = "mg/J"
  )
)
#TODO

Vos données sont maintenant correctement encodées.

Visualisation des données

Vous devez toujours débuter une analyse par visualiser vos données que ce soit sous la forme d’un tableau de données ou encore de graphiques.

Réalisation d’un tableau résumé des données

Pour réaliser un tableau résumé des données, vous avez à disposition le snippet .hmanova1desc qui renvoie les instructions suivantes :

DF %>.%
  group_by(., XFACTOR) %>.%
  summarise(., mean = mean(YNUM), sd = sd(YNUM), count = sum(!is.na(YNUM)))
toothgrowth <- read("ToothGrowth", package = "datasets", lang = "fr")
toothgrowth$dose <- as.ordered(toothgrowth$dose)
toothgrowth <- labelise(toothgrowth, self = FALSE,
  label = list(
    len = "Longueur des dents",
    supp = "Supplémentation",
    dose = "Dose"
  ),
  units = list(
    len = "mm",
    supp = NA,
    dose = "mg/J"
  )
)
toothgrowth %>.%
  group_by(., dose) %>.%
  summarise(., mean = mean(len), sd = sd(len), count = sum(!is.na(len)))
#TODO

Réalisation de graphiques

Vous trouverez dans la littérature plusieurs graphiques qui permettent de visualiser vos lorsque vous souhaitez réaliser une anova par la suite.

  • Réalisez une boite de dispersion. Le snippet est .cbbox et renvoie les instructions suivantes :
chart(data = DF, YNUM ~ XFACTOR) +
  geom_boxplot()
chart(data = toothgrowth, len ~ dose) +
  geom_boxplot()
#TODO

L’ANOVA a pour objectif de comparer des moyennes entre elles alors que la boîte de dispersion va vous représenter les 5 nombres qui sont des descripteurs non-paramétriques.

  • Réalisez un graphique en violon. Le snippet est .cbviolin et renvoit-e les instructions suivantes :
chart(data = DF, YNUM ~ XFACTOR) +
  geom_violin()
chart(data = toothgrowth, len ~ dose) +
  geom_violin() 
#TODO
  • Réalisez un graphique en violon en y ajoutant la moyenne de chaque groupe. Il n’existe pas de snippet réalisant directement ce graphique. Cependant, vous êtes capable de le réaliser par vous. Il s’agit d’une suite d’instructions.
chart(data = DF, YNUM ~ XFACTOR) + # Je réalise un graphique 
  geom_violin() + # Je réalise un graphique en violon
  geom_jitter(width = 0.05, alpha = 0.5) + # J'ajoute les points sur le graphique 
  stat_summary(geom = "point", fun.y = "mean", color = "red", size = 4) # J'ajoute la moyenne en rouge
chart(data = toothgrowth, len ~ dose) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.5) +
  stat_summary(geom = "point", fun.y = "mean", color = "red", size = 4)
#TODO

Maintenant que vous avez pris connaissance de vos données, vous pouvez réaliser votre test d’hypothèse en ayant pris connaissance de vos données.

Vérification des conditions d’applications

Test de Bartlett

Attention, avant de réaliser une analyse de variance vous devez vérifier si la variance de vos différents groupes est homogène au sein de la variable facteur que vous étudiez.

Le test de Bartlett est un outil mis à votre disposition. Le snippet est .hvbartlett et renvoie les instructions suivantes :

bartlett.test(data = DF, YNUM ~ XFACTOR)
bartlett.test(data = toothgrowth, len ~ dose)
#TODO

Distribution normale des résidus

Vous devez également vérifier sur vos résidus suivent une distribution normale.

  • Calculez vos résidus

Le snippet est .dvresid et renvoie les instructions suivantes :

DF %>.%
  mutate(., VAR.res = VAR - ave(VAR, FACTOR)) -> DF2
toothgrowth %>.%
  mutate(., len.res = len - ave(len, dose)) -> toothgrowth
#TODO
  • Réalisez un graphique quantile-quantile des résidus

Le snippet est .cuqqnorm et renvoie les instructions suivantes :

car::qqPlot(DF[["XNUM"]], distribution = "norm",
  envelope = 0.95, col = "Black", ylab = "XNUM")
toothgrowth <- read("ToothGrowth", package = "datasets", lang = "fr")
toothgrowth$dose <- as.ordered(toothgrowth$dose)
toothgrowth %>.%
  mutate(., len.res = len - ave(len, dose)) -> toothgrowth
car::qqPlot(toothgrowth[["len.res"]], distribution = "norm",
  envelope = 0.95, col = "Black", ylab = "len.res")
#TODO

Vos résidus suivent une distribution normale et la variance de vos trois groupes est homogènes. Vous pouvez maintenant réalisez votre ANOVA.

ANOVA

Réalisez l’analyse de variances. Le snippet est .hmanova1 et renvoie les instructions suivantes :

anova(anova. <- lm(data = DF, YNUM ~ XFACTOR))
anova(anova. <- lm(data = toothgrowth, len ~ dose))
#TODO

Une fois votre ANOVA réalisée, vous avez à votre disposition un snippet qui vous permet de vérifier la distribution normale de vos résidus. Cette méthode est plus rapide que de calculer vos résidus et d’en faire un graphique quantile-quantile. Cependant, il nécéssite de d’abord réaliser votre ANOVA et donc d’anticiper sur les conditions d’application. Néanmoins le test étant relativement robuste à de petites variations par rapport à la distribution normale, nous ne serons pas excessivement strict sur cette condition. On peut la vérifier après avoir réalisé l’anova.

#plot(anova., which = 2)
anova. %>.%
  chart(broom::augment(.), aes(sample = .std.resid)) +
  geom_qq() +
  #geom_qq_line(colour = "darkgray") +
  labs(x = "Theoretical quantiles", y = "Standardized residuals") +
  ggtitle("Normal Q-Q")
toothgrowth <- read("ToothGrowth", package = "datasets", lang = "fr")
toothgrowth$dose <- as.ordered(toothgrowth$dose)
anova. <- lm(data = toothgrowth, len ~ dose)

Vous observez que vous ne devez changer aucune instruction. Les snippets ont vraiment été conçu afin de simplifier vos analyses dans R.

#plot(anova., which = 2)
anova. %>.%
  chart(broom::augment(.), aes(sample = .std.resid)) +
  geom_qq() +
  #geom_qq_line(colour = "darkgray") +
  labs(x = "Theoretical quantiles", y = "Standardized residuals") +
  ggtitle("Normal Q-Q")
#TODO

Analyse Post Hoc

Votre ANOVA indique qu’il y a au moins une des moyennes qui est significativement différente des autres. Afin de connaitre le niveau dont la moyenne (ou les niveaux) qui est significativement différent des autres moyennes, vous pouvez réaliser une analyse complémentaire de l’ANOVA, une analyse post hoc dont le snippet est le suivant :

summary(anovaComp. <- confint(multcomp::glht(anova.,
  linfct = multcomp::mcp(XFACTOR = "Tukey")))) # Add a second factor if you want
.oma <- par(oma = c(0, 5.1, 0, 0)); plot(anovaComp.); par(.oma); rm(.oma)
summary(anovaComp. <- confint(multcomp::glht(anova.,
  linfct = multcomp::mcp(dose = "Tukey")))) # Add a second factor if you want
.oma <- par(oma = c(0, 5.1, 0, 0)); plot(anovaComp.); par(.oma); rm(.oma)
#TODO

Ayez une véritable réfléxion sur l’ensemble de l’analyse que vous venez de réaliser. Toutes les étapes de cette analyse sont importantes.

Conclusion

Bravo! Vous venez de terminer votre séance d’exercices dans un tutoriel “learnr”.

Laissez nous vos impressions sur cet outil pédagogique ou expérimentez encore dans la zone ci-dessous. Rappelez-vous que pour placer un commentaire dans une zone de code R, vous devez utilisez un dièse (#) devant vos phrases.

# Ajout de commentaires 
# ...
# Not yet...

Traitement des données I

Guyliann Engels & Philippe Grosjean