Gagnez du temps par l'automatisation de vos tests de permutation
Comme je l’expliquais dans un précédent article, les tests de permutation sont des approches robustes, basées sur l’aléatoire et le ré-échantilonnage, qui permettent de réaliser des tests statistiques paramétriques, sans que la validité des résultats ne reposent sur des distributions théoriques !
En pratique, les tests de permutations permettent, par exemple, de réaliser des régressions linéaires (simples ou multiples) sans avoir à vérifier les hypothèses de normalité, et d’homogénéité des
résidus. C’est tentant, non ?
Ces approches sont donc intéressantes, lorsque l’on est pressé par le temps et que l’on doit réaliser, de façon répétée, les mêmes modèles de régression, avec différentes variables “réponse”, par exemple. Mais pour aller encore plus vite, il peut être intéressant de réaliser une automatisation ces tests de permutation.
Dans ce post, je vais donc vous montrer comment automatiser la réalisation de régressions linéaires simples, à l’aide de tests de permutations. Pour cela, je vais utiliser le package `purrr` du package `tidyverse`, et le jeu de données `cystfibr` du package `ISwR`. Le package `ggplot2` (qui appartient aussi au package `tidyverse`) sera également utile, ainsi que le package `lmPerm` qui permet d’utiliser les tests de permutation.
Chargement des packages
library(tidyverse)
library(lmPerm)
library(ISwR) Les données
Les données s’intéressent aux fonctions pulmonaires de 25 patients atteints de fibrose cystique ; 10 variables sont mesurées.
data(cystfibr)
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 Nous allons automatiser la réalisation de régressions linéaires simples entre chaque variable de ce jeu de données et la variable pemax qui sera, dans tous les cas, la variable explicative. Pour plus de détail sur la signification des variables, vous pouvez consulter les informations du jeu de données en utilisant la commande suivante
?cystfibr Evaluation des relations linéaires
Les tests de permutations peremettent de s’affranchir de l’évaluation des hypothèses de normalité et d’homogénité des résidus. En revanche, il est toujours nécessaire de vérifier que la forme de la relation entre la variable réponse et la variable explicative est globalement linéaire.
Pour cela, on peut utiliser deux méthodes.
Passage en format long et utilisation de facet_wrap
Ici, il n’y a que 9 relations à évaluer, ce n’est pas beaucoup. On peut donc passer les données en format long à l’aide de la fonction gather, et utiliser la fonction facet_wrap de ggplot2.
cystfibr_long <- cystfibr %>%
gather(var, value, -pemax) Dans ce format long, les données des variables age à tlc sont placées les unes sous les autre dans une même colonne nommée var. Leur valeurs sont en vis à vis dans la colonne value. Et la variable pemax est répétée.
cystfibr_long[1:30,]
## pemax var value
## 1 95 age 7
## 2 85 age 7
## 3 100 age 8
## 4 85 age 8
## 5 95 age 8
## 6 80 age 9
## 7 65 age 11
## 8 110 age 12
## 9 70 age 12
## 10 95 age 13
## 11 110 age 13
## 12 90 age 14
## 13 100 age 14
## 14 80 age 15
## 15 134 age 16
## 16 134 age 17
## 17 165 age 17
## 18 120 age 17
## 19 130 age 17
## 20 85 age 19
## 21 85 age 19
## 22 160 age 20
## 23 165 age 23
## 24 95 age 23
## 25 195 age 23
## 26 95 sex 0
## 27 85 sex 1
## 28 100 sex 0
## 29 85 sex 1
## 30 95 sex 0 ggplot(cystfibr_long, aes(x=pemax, y=value))+
geom_point()+
geom_smooth( method=loess, se=FALSE)+
geom_smooth(method=lm, colour="red", fill="red", alpha=0.25) +
facet_wrap(~var, scales="free") +
theme_classic () 
Au final, la forme linéaire me semble acceptable pour toutes les relations. Celle relative à la variable `sex` n’a pas vraiment de sens, je vais donc la retirer du jeu de donnés qui sera fourni aux tests de permutation automatisés. Il faudrait faire de même avec les variables qui ne présentent pas de relation de forme linéaire avec la variable prédictive employée.
Réalisation automatisée des plots
Comme expliqué dans cet article, il est également possible de réaliser les plots un par un de façon automatique, et de les enregistrer dans le working directory (ou ailleurs). Pour cela, on crée notre propre fonction de plot, et on l’utilise dans une boucle.
plot_scatt <- function(df, x, y)
{
ggplot(df, aes_string(x=x, y=y))+
geom_point()+
geom_smooth( method=loess, se=FALSE)+
geom_smooth(method=lm, colour="red", fill="red",
alpha=0.25) +
theme_classic ()
}
for(i in 1 : (ncol(cystfibr)-1)) {
jpeg(paste(names(cystfibr)[i], "jpeg", sep = "."), width = 15, height =12,
units="cm", quality=75, res=300)
p <- plot_scatt(cystfibr,names(cystfibr)[i], "pemax")
print(p)
dev.off()
} Automatisation des tests de permutation
L’automatisation des tests de permutation se fait sur un jeu de données au format long.
cystfibr2_long <- cystfibr %>%
select(-sex) %>%
gather(var, value, -pemax) L’automatisation des tests de permutation consiste à:
- transformer le jeu de données long qui est un data.frame en une list avec la fonction `split`.
- utiliser la fonction `map` du package `purrr` pour appliquer le test de permutation sur chaque élement de la list.
- utiliser la fonction `summary`, toujours à l’aide de la **fonction `map`**, pour obtenir l’ensemble des résultats de chaque régression.
set.seed(1234)
res_complet <- cystfibr2_long %>%
split(.$var) %>%
map(~lmp(pemax~value, data=.)) %>%
map(summary) Les résultats sont alors stockés dans l’objet “res_complet”.
En ajoutant une fonction personnalisée aux lignes précédentes, il est alors possible de mettre en forme les résultats, de façon plus synthétique, comme ceci :
set.seed(1234)
res_resume <- cystfibr2_long %>%
split(.$var) %>%
map(~lmp(value~pemax, data=.)) %>%
map(summary) %>%
map(function(x) {
out <- as.data.frame(x$coefficients[2,, drop = FALSE])}) %>%
bind_rows(.id = "variable")
res_resume
## variable Estimate Iter Pr(Prob)
## 1 age 0.09281830 5000 0.00060000
## 2 bmp 0.08240561 295 0.25423729
## 3 fev1 0.15182405 3334 0.02939412
## 4 frc -0.54549981 777 0.11454311
## 5 height 0.38529940 5000 0.00040000
## 6 rv -0.81175762 630 0.13809524
## 7 tlc -0.09216387 90 0.53333333
## 8 weight 0.34002200 5000 0.00000000 Remarque : Pour mettre au point ce code, je me suis inspirée d’un post de Brian S. Yandell.
Voilà, j’espère que cet article vous sera utile, et qu’il vous permettra de gagner du temps lorsque vous aurez à faire des analyses répétées.
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édits photos : geralt
