Intervenant.e.s

Rédaction

Clément Poupelin, clementjc.poupelin@gmail.com

Relecture

Rappels sur les tests multiples

Setup

Show the code

# Données
library(dplyr)        # manipulation des données

# Plots
## ggplot
library(ggplot2)
library(gridExtra)

Simulation des p-values
Fonction qui calcul les résultats

Show the code

simulate_pvalues <- function(n, m, m0, mu1) {
  p_values <- numeric(m)
  
  for (i in 1:m0) {
    sample <- rnorm(n, mean = 0, sd = 1)  
    p_values[i] <- t.test(sample)$p.value
  }
  
  for (i in (m0+1):m) {
    sample <- rnorm(n, mean = mu1, sd = 1)  
    p_values[i] <- t.test(sample)$p.value
  }
  
  return(p_values)
}

Show the code

results_fun <- function(p_values, m0, alpha = 0.05) {
  positives <- p_values < alpha
  bonferroni_threshold <- alpha / length(p_values)
  bh_thresholds <- p.adjust(p_values, method = "BH")
  
  bonferroni_positives <- p_values < bonferroni_threshold
  bh_positives <- bh_thresholds < alpha
  
  count_results <- function(positives) {
    true_positives <- sum(positives[(m0 + 1):length(p_values)])
    false_positives <- sum(positives[1:m0])
    proportion_fp <- false_positives / max(1, sum(positives))
    return(c(
      sum(positives),
      true_positives,
      false_positives,
      proportion_fp
    ))
  }
  
  list(
    "Sans correction" = count_results(positives),
    "Bonferroni" = count_results(bonferroni_positives),
    "Benjamini-Hochberg" = count_results(bh_positives)
  )
}

Show the code

set.seed(140400)

Données

Dans cette exercice nous allons générer des p-values en différentes quantités afin d’illustrer l’importance des techniques de correction pour des tests mutliples.

On va alors utiliser la fonction simulate_pvalues qui prend en entrée $n$ , $m$ , $m_{0} \leq m$ et $μ_{1}$ , et qui :

simule $m_{0}$ échantillons contenant chacun $n$ réalisations d’une loi Normale centrée réduite, et $m - m_{0}$ échantillons de taille $n$ d’une loi normale d’espérance $μ_{1}$ et de variance $1$
effectue pour chacun des $m$ échantillons un test de Student de nullité de la moyenne
retourne les $m$ p-values associées à ces tests

Dans la suite, on considère que $n = 100$ et $m = 1000$ . Pour différentes situations, on va alors relever le nombre de positifs, de vrais-positifs, de faux-positifs et la proportion de faux-positifs lorsqu’on applique chacun des $m$ tests précédents au niveau $α = 0.05$ sans correction, lorsqu’on applique une procédure de Bonferroni associée à FWER $\leq 0.05$ et lorsqu’on applique une procédure de Benjamini-Hochberg associée à FDR $\leq 0.05$ .

Il y a donc pour chaque situation, 3 méthodes et 4 scores à calculer par méthode

Les différents scénarios testés seront :

Pas de positifs : $m_{0}$ = 1000
Peu de positifs, facilement identifiables : $m_{0}$ = 950, $μ_{1}$ = 1
Peu de positifs, difficilement identifiables : $m_{0}$ = 950, $μ_{1}$ = 0.3
Pas mal de positifs, facilement identifiables : $m_{0}$ = 800, $μ_{1}$ = 1
Pas mal de positifs, difficilement identifiables : $m_{0}$ = 800, $μ_{1}$ = 0.3
Beaucoup de positifs, facilement identifiables : $m_{0}$ = 200, $μ_{1}$ = 1
Beaucoup de positifs, difficilement identifiables : $m_{0}$ = 200, $μ_{1}$ = 0.3

Show the code

scenarios <- list(
  c(1000, 0),
  c(950, 1),
  c(950, 0.3),
  c(800, 1),
  c(800, 0.3),
  c(200, 1),
  c(200, 0.3)
)

Show the code

results <- list()
n <- 100
m <- 1000

for (i in 1:length(scenarios)) {
  m0 <- scenarios[[i]][1]
  mu1 <- scenarios[[i]][2]
  p_values <- simulate_pvalues(n, m, m0, mu1)
  results[[paste("m0 =", m0, "mu1 =", mu1)]] <- results_fun(p_values, m0)
}

# Création du tableau des résultats
results_matrix <- do.call(rbind, lapply(names(results), function(scenario) {
  cbind(Scenario = scenario, do.call(rbind, results[[scenario]]))
}))
results_df <- as.data.frame(results_matrix)
names(results_df) <- c("Scenario", "Total Positifs", "Vrais Positifs", "Faux Positifs", "Proportion Faux Positifs")

Analyse des scénarios

Ici on a pas de positifs : $m_{0}$ = 1000.

Show the code

results_df[1:3,] %>% DT::datatable()

Show entries

Search:

	Scenario	Total Positifs	Faux Positifs	Proportion Faux Positifs
Sans.correction	m0 = 1000 mu1 = 0	57	57	1
Bonferroni	m0 = 1000 mu1 = 0	0	0	0
Benjamini.Hochberg	m0 = 1000 mu1 = 0	0	0	0

Showing 1 to 3 of 3 entries

Previous1Next

Ici on a peu de positifs, facilement identifiables : $m_{0}$ = 950, $μ_{1}$ = 1

Show the code

results_df[4:6,] %>% DT::datatable()

Show entries

Search:

	Scenario	Total Positifs	Vrais Positifs	Faux Positifs	Proportion Faux Positifs
Sans.correction.1	m0 = 950 mu1 = 1	91	50	41	0.450549450549451
Bonferroni.1	m0 = 950 mu1 = 1	50	50	0	0
Benjamini.Hochberg.1	m0 = 950 mu1 = 1	51	50	1	0.0196078431372549

Showing 1 to 3 of 3 entries

Previous1Next

Ici on a peu de positifs, difficilement identifiables : $m_{0}$ = 950, $μ_{1}$ = 0.3

Show the code

results_df[7:9,] %>% DT::datatable()

Show entries

Search:

	Scenario	Total Positifs	Vrais Positifs	Faux Positifs	Proportion Faux Positifs
Sans.correction.2	m0 = 950 mu1 = 0.3	95	42	53	0.557894736842105
Bonferroni.2	m0 = 950 mu1 = 0.3	3	3	0	0
Benjamini.Hochberg.2	m0 = 950 mu1 = 0.3	8	8	0	0

Showing 1 to 3 of 3 entries

Previous1Next

Ici on a pas mal de positifs, facilement identifiables : $m_{0}$ = 800, $μ_{1}$ = 1

Show the code

results_df[10:12,] %>% DT::datatable()

Show entries

Search:

	Scenario	Total Positifs	Vrais Positifs	Faux Positifs	Proportion Faux Positifs
Sans.correction.3	m0 = 800 mu1 = 1	240	200	40	0.166666666666667
Bonferroni.3	m0 = 800 mu1 = 1	200	200	0	0
Benjamini.Hochberg.3	m0 = 800 mu1 = 1	207	200	7	0.0338164251207729

Showing 1 to 3 of 3 entries

Previous1Next

Ici on a pas mal de positifs, difficilement identifiables : $m_{0}$ = 800, $μ_{1}$ = 0.3

Show the code

results_df[13:15,] %>% DT::datatable()

Show entries

Search:

	Scenario	Total Positifs	Vrais Positifs	Faux Positifs	Proportion Faux Positifs
Sans.correction.4	m0 = 800 mu1 = 0.3	208	166	42	0.201923076923077
Bonferroni.4	m0 = 800 mu1 = 0.3	25	25	0	0
Benjamini.Hochberg.4	m0 = 800 mu1 = 0.3	128	123	5	0.0390625

Showing 1 to 3 of 3 entries

Previous1Next

Ici on a beaucoup de positifs, facilement identifiables : $m_{0}$ = 200, $μ_{1}$ = 1

Show the code

results_df[16:18,] %>% DT::datatable()

Show entries

Search:

	Scenario	Total Positifs	Vrais Positifs	Faux Positifs	Proportion Faux Positifs
Sans.correction.5	m0 = 200 mu1 = 1	810	800	10	0.0123456790123457
Bonferroni.5	m0 = 200 mu1 = 1	800	800	0	0
Benjamini.Hochberg.5	m0 = 200 mu1 = 1	809	800	9	0.0111248454882571

Showing 1 to 3 of 3 entries

Previous1Next

Ici on a beaucoup de positifs, difficilement identifiables : $m_{0}$ = 200, $μ_{1}$ = 0.3

Show the code

results_df[19:21,] %>% DT::datatable()

Show entries

Search:

	Scenario	Total Positifs	Vrais Positifs	Faux Positifs	Proportion Faux Positifs
Sans.correction.6	m0 = 200 mu1 = 0.3	682	668	14	0.0205278592375367
Bonferroni.6	m0 = 200 mu1 = 0.3	91	91	0	0
Benjamini.Hochberg.6	m0 = 200 mu1 = 0.3	631	622	9	0.0142630744849445

Showing 1 to 3 of 3 entries

Previous1Next

Résultats

Les avantages et inconvénients des méthodes sont :

Sans correction : Forte sensibilité mais haut risque de faux positifs
Bonferroni : Réduit considérablement les faux positifs mais manque de puissance
Benjamini-Hochberg : Bon compromis entre puissance et contrôle des faux positifs

Conclusion

Session info

Show the code

sessioninfo::session_info(pkgs = "attached")

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.2 (2024-10-31)
 os       Ubuntu 24.04.1 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2025-02-27
 pandoc   3.2 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package   * version date (UTC) lib source
 dplyr     * 1.1.4   2023-11-17 [1] CRAN (R 4.4.2)
 ggplot2   * 3.5.1   2024-04-23 [1] CRAN (R 4.4.2)
 gridExtra * 2.3     2017-09-09 [1] CRAN (R 4.4.2)

 [1] /home/clement/R/x86_64-pc-linux-gnu-library/4.4
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────