Réalisez de multiples tests en quelques lignes avec une fonction map()

Cette semaine, au cours d’une séance de coaching “libre”, on m’a demandé comment faire de multiples tests de wilcoxon de façon automatisés (il y en avait plus de 250 à faire), et comment récupérer la p-value, afin ensuite de faire un sous jeu de données avec seulement les variables qui seraient significatives.

J’ai trouvé cette question vraiment intéressante, parce que c’est vrai que ça peut être super utile de savoir le faire ! Et aussi parce que j’y ai vu une bonne occasion de me (re)pencher sur les fonctions `map()` du packages `purrr` dont je parlais dans cet article celles du package`broom`.

Prérequis

 Pour tirer pleinement profit cet article il est nécessaire de connaitre l’utilisation du package dplyr et du symbole pipe `%>%` qui permet d’enchaîner les instructions.

Le jeu de données utilisé

Pour illustrer cet article, je vais reprendre les données `cystfibr` que j’avais utilisées dans l’article “Gagnez du temps en utilisant des tests de permutation et en les automatisant !” dans lequel je parlais déjà d’automatisation. Ce jeu de données `cystfibr` est un data frame de 25 lignes et 10 variables. Il est relatif aux fonctions pulmonaires de patients atteint de fibrose cystique. `cystfibr` est contenu dans le package `IswR` :
library(ISwR)
head(cystfibr)

    ##   age sex height weight bmp fev1  rv frc tlc pemax
    ## 1   7   0    109   13.1  68   32 258 183 137    95
    ## 2   7   1    112   12.9  65   19 449 245 134    85
    ## 3   8   0    124   14.1  64   22 441 268 147   100
    ## 4   8   1    125   16.2  67   41 234 146 124    85
    ## 5   8   0    127   21.5  93   52 202 131 104    95
    ## 6   9   0    130   17.5  68   44 308 155 118    80 

Le test de Wilcoxon

Le test de Wilcoxon est un test non paramétrique (basé sur le calcul d’une somme de rangs) qui permet de comparer les médianes de deux groupes. Avec le jeu de données `cystfibr`, il s’agit, par exemple, de tester si les médianes des pressions expiratoires maximales (variables pemax) sont différents entre les hommes (n=14) et les femmes (n=11).

Le test de wilcoxon se réalise avec la commande suivante :

wilcox.test(pemax~sex, data=cystfibr)

    ## 
    ##  Wilcoxon rank sum test with continuity correction
    ## 
    ## data:  pemax by sex
    ## W = 100, p-value = 0.2158
    ## alternative hypothesis: true location shift is not equal to 0 

La sortie comporte plusieurs informations :

  • un rappel des données employées
  • la statistique du test : W
  • la p-value du test
  • un rappel de l’hypothèse alternative (la différence des médianes est différente de 0)

Réalisez de multiples tests de façon automatisée

Il s’agit de réaliser un test de Wilcoxon sur chacune des neuf variables numériques, pour comparer les médianes relatives aux hommes et aux femmes. Autrement dit, il s’agit de réaliser neuf tests de Wilcoxon. Neuf tests de Wilcoxon, ce n’est pas la mer à boire. Mais dans la demande initiale qui m’a été faite il y en avait près de 250 !

Multiples tests : approche simple

L’approche simple consiste à :

  • supprimer la variable qui contient les deux groupes à comparer (ici sex)
  • utiliser la fonction ‘map()’ (du package ‘purrr’) pour réaliser un test de wilcoxon avec chacune des variables numérique et la variable sex, qui elle est contenue dans le jeu cystfibr (c’est pour cela que l’on indique ‘data=cystfibr’). A ce stade toutes les sorties de chaque test de wilcoxon sont stockées dans une liste.
  • utiliser la fonction ‘map_dbl()’ pour extraire la p.value de la sortie (voir plus haut). La fonction ’map_dbl()’ permet de concaténer l’ensemble des p-values dans un vecteur (nommé ici p_values)
library(tidyverse)
library(purrr)

pvalues <- cystfibr%>% 
      select(-sex) %>% 
      map(~wilcox.test(.x~sex, data=cystfibr)) %>% 
      map_dbl("p.value")


pvalues <- round(pvalues,3) # arrondi à 3 chiffres après la virgule 

Dans un second temps, on récupère le nom des variables numériques :

var <- names(cystfibr)[c(1,3:10)] 

Puis on peut stoker les nom de ces variables et leur p-value correspondante, dans un data frame :

pval_df <- data_frame(var,pvalues)
pval_df

    ## # A tibble: 9 x 2
    ##   var    pvalues
    ##       
    ## 1 age      0.441
    ## 2 height   0.261
    ## 3 weight   0.344
    ## 4 bmp      0.493
    ## 5 fev1     0.009
    ## 6 rv       0.218
    ## 7 frc      0.476
    ## 8 tlc      0.805
    ## 9 pemax    0.216 
Il est alors tout à fait possible d’ajuster les p-values pour prendre en compte cette réalisation de multiples tests et le problème de l’augmentation du risque alpha global. Ici, avec la méthode d’ajustement de Benjamini & Hochberg :
pval_df$pval_ajust <- p.adjust(pval_df$pvalues, method="BH")

pval_df

    ## # A tibble: 9 x 3
    ##   var    pvalues pval_ajust
    ##             
    ## 1 age      0.441     0.555 
    ## 2 height   0.261     0.555 
    ## 3 weight   0.344     0.555 
    ## 4 bmp      0.493     0.555 
    ## 5 fev1     0.009     0.0810
    ## 6 rv       0.218     0.555 
    ## 7 frc      0.476     0.555 
    ## 8 tlc      0.805     0.805 
    ## 9 pemax    0.216     0.555 

Multiples tests : approche plus complexe

L’approche décrite ci-dessous est plus complexe, mais elle vaut la peine car elle fait appelle à des fonctions très intéressantes des packages `tidyr`, comme la fonction `nest()`, et du package `broom`, comme la fonction `tidy()`. De plus, cette approche peut très facilement être modifiée pour, par exemple, stocker les sorties de modèles de régressions linéaires (les coefficients, leur erreur standard, les pvalues, le r² etc…) Il s’agit dans un premier temps de réaliser une version “long” du jeu de données, c’est à dire de le réduire à 3 colonnes :
  • une contenant la variable “sex”,
  • une contenant toutes les valeurs numériques,
  • une contenant le nom de la variable à laquelle la valeur numérique se rapporte.
cystfibr_long <- cystfibr %>% 
      gather(var, value, -sex) 

head(cystfibr_long)

    ##   sex var value
    ## 1   0 age     7
    ## 2   1 age     7
    ## 3   0 age     8
    ## 4   1 age     8
    ## 5   0 age     8
    ## 6   0 age     9 

Dans un deuxième temps, on crée une version “nested” de la version “long”. Il s’agit de découper le jeu de données par paquet de lignes correspondant à chacune des variables, et de stocker les données à proprement parlé, dans une variable appelée “data”. Ça à l’air un peu compliqué, mais on s ’y habitue…

cystfibr_nested <- cystfibr_long%>% 
      group_by(var) %>% 
      nest()  

head(cystfibr_nested )

    ## # A tibble: 6 x 2
    ##   var    data             
    ##                
    ## 1 age    
    ## 2 height 
    ## 3 weight 
    ## 4 bmp    
    ## 5 fev1   
    ## 6 rv 

Vous pouvez visualiser les données “nested” dans le viewer:

View(cystfibr_nested) 

En cliquant sur une case blanche de la colonne “data”, il est possible de visualiser le jeu de données qu’elle contient :

automatiser les tests statistiques avec le logiciel R

Ensuite, il s’agit de créer une fonction pour réaliser le test de wilcoxon de manière spécifique avec nos données, c’est à dire en spécifiant le nom des variables, tel qu’elles seront nommées dans le jeu de données (ici “value” et “sexe) :

wilcox_test <- function(df){
        wilcox.test(value~sex, data=df)
    } 

Puis de créer une nouvelle variable dans le data.frame “cystfibr_nest” qui contiendra la sortie du test de Wilcoxon créé à l’étape précédente :

cystfibr_nested <- cystfibr_nested %>% 
      mutate(wilcox_test=map(data,wilcox_test)) 

On utilise alors la fonction `tidy()` du package `broom` pour formater les sorties des tests de wilcoxon sous forme de data frame. Voici un exemple :

    wt <- wilcox.test(pemax~sex, data=cystfibr)
    wt

    ## 
    ##  Wilcoxon rank sum test with continuity correction
    ## 
    ## data:  pemax by sex
    ## W = 100, p-value = 0.2158
    ## alternative hypothesis: true location shift is not equal to 0

    broom::tidy(wt)

    ## 
    # A tibble: 1 x 4
    ##   statistic p.value method                                      alternative
    ##                                                        
    ## 1       100   0.216 Wilcoxon rank sum test with continuity cor~ two.sided 

Cette fonction `tidy()` est couplée à la fonction `map()` pour être appliquée à toutes les lignes (de cystfibr_nested), et à `mutate()` pour que le data frame issu de la fonction `tidy()` soit inclus dans une nouvelle variable (nommé tidy !) :

cystfibr_nested <- cystfibr_nested %>% 
      mutate(tidy=map(wilcox_test, broom::tidy)) 

Voici à quoi ressemble à présent le jeu de données “cystfibr_nested”

Et ici, le data frame stocké dans la variable tidy :

Il ne reste plus qu’à extraire les variables de ce data frame contenu dans la colonne “tidy” pour les ajouter directement aux colonnes du data frame “cystfibr_nested”. Ceci est possible avec la fonction `unnest()` :

 cystfibr_nested <- cystfibr_nested %>% 
      unnest(tidy, .drop=TRUE) 

L’argument `.drop=TRUE`, permet de supprimer la variable `tidy`.

cystfibr_nested

    ## # A tibble: 9 x 5
    ##   var    statistic p.value method                               alternative
    ##                                                   
    ## 1 age         91.5 0.441   Wilcoxon rank sum test with continu~ two.sided  
    ## 2 height      98   0.261   Wilcoxon rank sum test with continu~ two.sided  
    ## 3 weight      95   0.344   Wilcoxon rank sum test               two.sided  
    ## 4 bmp         90   0.493   Wilcoxon rank sum test with continu~ two.sided  
    ## 5 fev1       126.  0.00856 Wilcoxon rank sum test with continu~ two.sided  
    ## 6 rv          54   0.218   Wilcoxon rank sum test with continu~ two.sided  
    ## 7 frc         63.5 0.476   Wilcoxon rank sum test with continu~ two.sided  
    ## 8 tlc         72   0.805   Wilcoxon rank sum test with continu~ two.sided  
    ## 9 pemax      100   0.216   Wilcoxon rank sum test with continu~ two.sided 

Finalement, en enchaînant les commandes, on obtient :

wilcox_test <- function(df){
        wilcox.test(value~sex, data=df)
    } 
wilcox_res <- cystfibr %>% 
      gather(var, value, -sex) %>% 
      group_by(var) %>% 
      nest()  %>% 
      mutate(
        wilcox_test=map(data,wilcox_test),
        tidy=map(wilcox_test, broom::tidy)) %>% 
      unnest(tidy,.drop=TRUE) 

On peut là encore ajuster les p-values :

wilcox_res$pval_ajust <- p.adjust(wilcox_res$p.value, method="BH")

    wilcox_res

    ## # A tibble: 9 x 6
    ##   var    statistic p.value method                    alternative pval_ajust
    ##                                              
    ## 1 age         91.5 0.441   Wilcoxon rank sum test w~ two.sided       0.554 
    ## 2 height      98   0.261   Wilcoxon rank sum test w~ two.sided       0.554 
    ## 3 weight      95   0.344   Wilcoxon rank sum test    two.sided       0.554 
    ## 4 bmp         90   0.493   Wilcoxon rank sum test w~ two.sided       0.554 
    ## 5 fev1       126.  0.00856 Wilcoxon rank sum test w~ two.sided       0.0770
    ## 6 rv          54   0.218   Wilcoxon rank sum test w~ two.sided       0.554 
    ## 7 frc         63.5 0.476   Wilcoxon rank sum test w~ two.sided       0.554 
    ## 8 tlc         72   0.805   Wilcoxon rank sum test w~ two.sided       0.805 
    ## 9 pemax      100   0.216   Wilcoxon rank sum test w~ two.sided       0.554 

Puis finalement faire un subset du tableau de résultat en ne gardant que les p-values inférieures à un certain seuil, 0.1 par exemple :

 wilcox_res_sig <- wilcox_res %>% 
      select(var, p.value, pval_ajust) %>%  # on ne garde que ces 3 var
      filter(pval_ajust<0.1)

wilcox_res_sig

    ## # A tibble: 1 x 3
    ##   var   p.value pval_ajust
    ##            
    ## 1 fev1  0.00856     0.0770 

Ajuster de multiples modèles

Si l’on souhaite, par exemple, faire toutes les régressions linéaires simples entre la variable “pemax” et chacune des autres variables numériques :

mod_lm <- function(df){
        lm(pemax~value, data=df)
    }

rls_res <- cystfibr %>% 
      gather(var, value, -sex,-pemax) %>% 
      group_by(var) %>% 
      nest()  %>% 
      mutate(
        mod=map(data, mod_lm),
        glance=map(mod, broom::glance)) %>% 
      unnest(glance,.drop=TRUE) %>% 
      select(var, r.squared, p.value)%>% 
      select(var,p.value) %>% 
      mutate(pval_ajust=p.adjust(p.value, method="BH"))  

rls_res

    ## # A tibble: 8 x 3
    ##   var     p.value pval_ajust
    ##              
    ## 1 age    0.00111     0.00413
    ## 2 height 0.00155     0.00413
    ## 3 weight 0.000646    0.00413
    ## 4 bmp    0.270       0.308  
    ## 5 fev1   0.0228      0.0457 
    ## 6 rv     0.124       0.166  
    ## 7 frc    0.0380      0.0608 
    ## 8 tlc    0.385       0.385 

Conclusion

J’espère que cet article vous aura convaincu de l’intérêt de savoir utiliser les fonctions `map()` du package `purrr`, mais aussi de l’intérêt de les coupler aux fonctions du package `broom`. Pour aller plus loin, je vous conseille de visionner cet talk d’Hadley Wickham :

Et pour découvrir le package broom, vous pouvez consulter sa vignette, ici

Si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee 🙏

Crédit photos : alan9187

Une réponse

Laisser un commentaire

Votre adresse de messagerie ne sera pas publiée. Les champs obligatoires sont indiqués avec *