Tests multiples - procédure par bootstrap

Bonjour,

Selon la méthodo que l'on peut retrouver dans ce document: http://statistics.berkeley.edu/sites/default/files/tech-reports/633.pdf p.14-15 (cf box1, ici on présente le cas par permutation mais normalement celui par bootstrap est similaire), voici le programme que j'ai écris (et qui évidemment ne marche pas hein):

Code:

bdd.multtest = function(BDD) {

 # Les procédures d'ajustement par permutation, SD permutation, FDR permutation, bootstrap, SD bootstrap, FDR bootstrap sont implémentées pour cette fonction

 N = dim(BDD)[1] # nombre d'observations
 P = dim(BDD)[2] # nombre de variables sachant que la première est la variable classe
 biblio.Y = summary(BDD[,1]) # descriptif de la variable classe
 nbClass.Y = length(biblio.Y) # nombre de classe de la variable classe
 B = 3000 # fixé à 3000 rééchantillonnage par bootstrap
 part = 1 # on tire autant de fois qu'on a d'observations
 BDD = data.frame(Y = BDD[,1],scale(BDD[,2:P])) # création de la base de données centrée-réduite
 
 # Création de la matrice de résultats
 p.correct = matrix(0,P-1,7) # initialisation de la matrice de résultats finaux pour les 7 méthodes prévues
 for (p in 2:P) { # pour chaque variable continue on calcul le test de Student pour la variable classe versus les différentes variables continues
 T = as.data.frame(cbind(BDD[,1],BDD[,p]))
 names(T) = c("Y","X")
 p.correct[p-1,1] = t.test(X~Y,T,var.equal=TRUE)$p.value # test de Student
 }

 # On applique la méthode par bootstrap
 p.b = matrix(0,P-1,B) # initialisation de la matrice contenant les différentes p-valeurs calculées pour le test de Student sur les B sous-échantillons construit par bootstrap
 for (b in 1:B) { # application du bootstrap
 ind.cl1 = BDD[which(BDD[,1]==names(biblio.Y)[1]),] # on initialise la procédure en traitant le cas particulier pour la première classe de la variable classe
 X.samp = cbind(ind.cl1[sample(1:dim(ind.cl1)[1],floor(part*biblio.Y[1]),replace=TRUE),]) # boostrap sur la classe 1
 for (k in 2:nbClass.Y) { # pour les k>1 autres classes on procède automatiquement à la même procédure que pour la classe 1
 ind.clk = BDD[which(BDD[,1]==names(biblio.Y)[k]),]
 samp.add = ind.clk[sample(1:dim(ind.clk)[1],floor(part*biblio.Y[k]),replace=TRUE),]
 X.samp = rbind(X.samp,samp.add)
 }
 X.samp = as.data.frame(X.samp) # edition de la matrice contenant le sous-échantillon pour le bème bootstrap
 for (p in 2:P) { # pour les p variables continues on calcul le test de Student pour le bème sous-échantillon construit par bootstrap
 X.test = X.samp[,c(1,p)]
 names(X.test) = c("Y","X")
 p.b[p-1,b] = t.test(as.numeric(X)~Y,X.test,var.equal=TRUE,alternative = "two.sided")$p.value # test de Student
 }
 }
 for (p in 2:P) {p.correct[p-1,2] = sum(p.b[p-1,] <= p.correct[p-1,1])/B} # les B bootstrap ayant eu lieu, on calcul les p-valeurs corrigées par procédure bootstrap selon la théorie
 return (p.correct)
}

Que fais ce code? Dans un premier temps il calcule le test de Student pour les différents Xi afin de déterminer les p-valeurs bruts. Ensuite, pour B tirage avec remise il construit la matrice de données par bootstrap contenant tous les Xi (en effet il faut faire un bootstrap général pour conserver l'état des corrélations entre les tests). Pour chacune de ces matrices on calcul le test de student qu'on stocke dans p.b. Une fois les B tirages effectués, on applique l'agorithme qui consiste à déterminer le nombre de p-valeurs calculées inférieures ou égales à la p-valeur brut, ce cardinal étant diviser par B assez logiquement.

En PJ j'ai mis un exemple basique et pour lequel le programme marche.

Voici le résultat que j'obtiens:

            pBrut          pBoot        SolSAS
7.139178e-07  0.6233333    <0.0001
8.128818e-01  0.8196667         1
6.942499e-01  0.7023333      0.9940
4.086682e-04  0.5416667       0.047
6.654268e-04  0.5676667       0.053

La première colonne contient les p-valeurs bruts du test de Student pour les 5 variables continues (donc sur l'échantillon complet/d'origine) et en seconde colonne leur ajustement par bootstrap. La troisième colonne représente la solution par SAS (que je ne comprends pas puisqu'on lorsqu'on jète un oeil à la matrice qu'il conçoit par boostrap et qui contient pour chaque tirage le sous-échantillon associé, mais mélange les observations quelque soit la classe de Y. Par exemple, on se retrouve avec des X pour Y = classe 1 dans la classe 2 alors que pour moi la définition du boostrap est un tirage avec remise en respectant leur appartenance aux différentes classes).

Personnellement je suis perdu, alors il est vrai que si on diminune le paramètre "part" de la fonction et qui correspond au nombre de tirage que l'on fait pour chaque classe on converge petit à petit vers un résultat cohérent mais jamais proche de celui de SAS (notamment pour la 3ème p-valeur du tableau ci-dessus où ça varie autour de 0.7 dans mon code et jamais autour de 1 comme sous SAS). A noter que SAS procède avec ce paramètre initialisé identiquement de mon code. Ci-dessous le code SAS utilisé également pour faire gagner du temps à l'âme charitable qui se sent de me filer un coup de main:

Code:


proc import datafile = "TEST.txt" out = TEST; run; /* importation de la base de données en format txt */

proc multtest data=test bootstrap nsample=3000 outsamp=res; /* option procédure bootstrap pour la proc multtest */
   test mean(X1 X2 X3 X4 X5); /* test de student pour X1, X2, X3, X4, X5 */
   class Y; /* précision de la variable classe */
run;

proc export data = res outfile="VUE.txt"; run; /* export de la matrice contenant les différents sous-échantillons construient par bootstrap */


A noter que lorsque je récupère la matrice contenant les différents sous-échantillons construient par bootstrap et que je calcule les différents tests de Student, je ne retombe pas sur la même chose que SAS quand j'applique l'algorithme... Je suis paumé nanananère...

Une idée? Car je suis proche de la corde et du tabouret là, d'ailleurs ils font des promos à Alinéa en ce moment...


Cordialement,

Salut,

tu trouveras peut-être des réponses sur ce site que j'ai découvert presque par hasard :

https://lemakistatheux.wordpress.com/category/tests-multiples/

hihihi

Niaboc

Ah et tu te crois drôle en prime

D'ailleurs ton avis sur ce site dont je ne vois absolument pas du tout qui est l'auteur.

Mon avis sur le site?

les plus :
+ très très large panoplie de tests déjà écrits (et encore plus à venir).
+ s'adresse aussi bien à des initiés en statistiques qu'à des noobs.
+ bonne maîtrise du sujet.
+ bien écrit.
+ un exemple concret sur tous les tests.
+ globalement, c'est un site à connaître et qui va rester dans les favoris.

les moins :
- trop centré sur les tests, il manque tout ce qu'on peut utiliser énormément en stat : les parties analyses factorielles, régressions, etc. qui viendront plus tard peut-être?
- trop de tests... il y a 50% des tests cités qui ne serviront jamais dans un cadre professionnel "classique". Quelqu'un d'impatient ne passera pas son temps à chercher ce qu'il veut... Il faudrait peut-être classer les tests selon leur utilité?

Mais un très bonne impression quand même.
On se demande comment l'auteur a le temps d'écrire tout ça? Entre deux parties de Star Wars Galaxy peut-être?

Merci pour ta réponse

Ce site étant ma première expérience et c'est clair que j'ai vite appris à mes dépends que ça représente un travail à plein temps surtout quand on veut viser l'exhaustivité et qu'au final on se rend compte qu'on a facilement entre 150 et 200 outils de statistiques toute application voulue confondue... Alors du coup j'ai rédigé la partie élaboration d'un projet qui me sert également de planning, donc au long et très long terme il y aura également la partie imputation des données manquantes, calcul des pondérations, analyse exploratoire, analyse supervisée et non supervée. Le site sera également étendu dans un temps très très éloigné à la théorie de la mesure et les algorithmes stochastiques et de la bioinformatique ou du moins les outils qui lui sont utiles, ainsi qu'aux paradoxes probabilistes (paradoxe du cercle et du triangle, de la fin du monde, de la vie sur une autre planète etc) et des projets que j'ai en tête comme le PMU ou encore la typologie de population à Mayotte. J'avais même pensé un temps à faire un petit truc lié au fofo sur l'évolution du nombre de topics par mois avec les séries chronologiques ou modèles de Markov avant de me souvenir que je delete les posts qui ont deux mois d'ancienneté et n'ont pas de réponse.

D'ailleurs dans la partie brouillon (lecture à risque puisque non relu) il y a l'ACM, l'ACP, la régression logistique et la version PLS, les K plus proches voisins pondérés (approche discriminante), l'analyse discriminante de Fisher linéaire et quadratique, les arbres de décision et les forêts aléatoire. Je rêve encore de pouvoir mettre les pages sur les SVM et les réseaux de neurones mais je n'ai encore jamais réussi à faire tourner un exemple de A à Z qui puisse confirmer la théorie...

Pour la petite histoire, et totalement inattendue, la page qui est la plus regardée est celle du test de variance de Fisher-Snedecor avec un total de 8 715 vues... suivi du test de Shapiro-Wilk avec 1 532 vues. Donc pour le moment je finalise les pages en fonction du nombre de vue, ce qui rejoins l'une de tes remarques. C'est clair que les tests associés aux séries chronologiques auront à tout casser une à deux vues par année pour la plupart.

Niveau ergonomie il faut que je prenne le temps de tester les différents thèmes gratuit et payant de wordpress mais plus ça avance plus je me rend compte que ça va vite devenir le bordel comme présentation, surtout qu'à la base les différents articles appartiennent à des catégories que j'ai créé mais qui n'apparaissent pas avec le thème de présentation que j'utilise en ce moment.

Pas mal la référence sur Star War Galaxies (mes plus belles nuits blanches quand j'étais un étudiant geek, évidemment la honte m'empêche de parler de celles sur World of Warcraft, toujours eu un doute à inscrire cette période de ma vie dans mon CV)   T'as chopé ça ou? Pitié pas jeuxvideo.com, la honte quand je pense au défouloire de certains de mes posts sur fond d'écriture sms kikoolol  

Hello,
Je viens de decouvrir le site. Super initiative, et super clair! Bon courage pour continuer ce gigantesque chantier (sans fin?). Je le rajouterai dans la liste des sites a consulter pour mes etudiants!

Ayana

joyeux_lapin13 a écrit:T'as chopé ça ou? Pitié pas jeuxvideo.com, la honte quand je pense au défouloire de certains de mes posts sur fond d'écriture sms kikoolol  

Merci Ayana

Oh putain Niaboc Quand le passé te rattrape...