Test de permutation

Jaromil Frossard
octobre 2015

Introduction

Méthode :

  1. Calculer une statistique (\( F_Y \)).
  2. Calculer cette statistique sous les données permutées.
  3. Répéter 2. afin d'obtenir une distribution \( \hat{G}_{F_Y} \) sous \( H_0 \) de la statistique.
  4. Comparer \( F_Y \) avec \( \hat{G}_{F_Y} \). La proportion de valeur de \( \hat{G}_{F_Y} \) plus extrême que \( F_Y \) est interprêtée comme une p-valeur.

Intuition : La permutation détruit la relation entre la VD et la VI; on produit donc des données sous \( H_0 \).

Introduction

Posulats et intérêts

Sous \( H_0 \) les observations doivent être échangeables. La distribution multivariée de \( Y \) est identique après permutation.

Le test est exact : sous \( H_0 \), on rejette \( H_0 \) exactement à un niveau \( \alpha \). La p-valeur du test sous \( H_0 \) suit une loi \( \mathcal{U}(0,1) \).

Les p-valeurs triées comparées à leur rang forment une droite.

Limitation

Dans un modèle d'ANOVA factorielle, les différents facteurs empêchent la création de la bonne distribution sous \( H_0 \).

ex.:

\( Y\sim \textrm{sexe}*\textrm{traitement} \)

Si on permute \( Y \), on obtient un distribution sous

\( H_0 : \beta_{ \textrm{sexe}}=0 \) et \( \beta_{ \textrm{traitement}}=0 \) et \( \beta_{ \textrm{sexe * traitement}}=0 \).

On ne peut donc pas tester indépendamment les variables sexe et traitement. Ce même problème apparaît aussi pour l'ANCOVA ou la RLM.

Formalisation du problème

\[ Y=X\beta +\epsilon = X_1\beta_1+X_2\beta_2+\epsilon \]

\( Y \) : VD (dim: \( n \times 1 \)),

\( X_1 \) : VI de nuissance (dim : \( n \times q \)),

\( X_2 \) : VI d'intérêt (dim : \( n \times (p-q) \)),

\( \epsilon \) erreur (dim: \( n \times 1 \)).

Statistique de test :

\[ F=\frac{RSS_{\textrm{petit}}-RSS_{\textrm{grand}}}{RSS_{\textrm{grand}}}\frac{n-p}{p-q} \]

Matrices "hat" et de résidus

\( H_X \): matrice qui produit les prédictions d'une régression en fonction de \( X \). \( H_X= X(X^\top X)^{-1}X^\top \).

\( R_X \): matrice qui produit les résidus d'une régression en fonction de \( X \). \( R_X= I-X(X^\top X)^{-1}X^\top \)

exemple : \( R_X X=0 \), ou \( (H_X+R_X)Y=Y \)

Pour tester \( H_0 : \beta_2=0 \)

\[ F=\frac{(R_{X_1}Y)^\top (R_{X_1}Y) -(R_{X}Y)^\top (R_{X}Y)}{(R_{X}Y)^\top (R_{X}Y)}\frac{n-p}{p-q} \]

Matrice de permutation

Une permutation peut s'écrire sous forme matricielle en considérant \( P \), qui est une matrice est composée de 1 ou de 0. La somme de chaque ligne et de chaque colonne vaut 1.

exemple: \[ P=\begin{pmatrix}1 &0&0 \\ 0&0&1\\0&1&0 \end{pmatrix}, \quad \textrm{ou} \quad P\begin{pmatrix}y_1 \\ y_2 \\ y_3 \end{pmatrix} =\begin{pmatrix}y_1 \\ y_3 \\ y_2 \end{pmatrix} \]

On permute le vecteur \( Y \) en calculant \( PY \) et la matrice \( M \) avec \( PMP^{\top} \).

Contrôle de la variable de nuisance

Paramétrique \[ Y= X_1\beta_1+X_2\beta_2+\epsilon \textrm{, ou } \epsilon\sim \mathcal{N}(0,\sigma^2I) \]

Kennedy \[ PR_{1}Y= R_{1}X_2\beta_2+R_{1}\epsilon \]

Freedman-Lane \[ (H_1+PR_1)Y= X_1\beta_1+X_2\beta_2+\epsilon \]

Huh-Jhun \[ PVR_{1}Y= VR_{1}X_2\beta_2+VR_{1}\epsilon \]

Simulation Monte-Carlo : le modèle

\[ Y=X\beta +\epsilon \]

\( Y \) VD (dim: \( n \times 1 \)),

\( X \) VI (dim: \( n \times 4 \)) composé de intercept, \( \mathcal{N}(0,1) \), \( \mathcal{N}(0,1) \) et \( \mathcal{U}(0,1) \),

\( \epsilon \sim \mathcal{N}(0,1) \) erreur (dim: \( n \times 1 \)).

\( \beta= (3 \quad 0 \quad 4 \quad -1.5)^\top \), \( n=30 \), \( 2000 \) réplications

\( \beta_2 \) est le paramètre sous \( H_0 \).

Simulation Monte-Carlo

plot of chunk unnamed-chunk-1

Simulation Monte-Carlo

plot of chunk unnamed-chunk-2

Statistique robuste

Permutations de type Freedman-Lane

\[ PW^{1/2}P^\top \mathbf{(H^W_1+PR^W_1)Y}= \\ PW^{1/2}P^\top \mathbf{X_1\beta_1}+PW^{1/2}P^\top \mathbf{X_2\beta_2}+PW^{1/2}P^\top \mathbf{\epsilon} \]

ou \( W \) sont des poids estimés sur le petit modèle.

\( H^W_1 \) et \( R^W_1 \) sont respectivement les matrices “hat” et “de résidus” associées à une estimation robuste du petit modèle, \( Y \sim X_1 \).

La pré-multiplication par \( PW^{1/2}P \) permet d'obtenir des RSS pondérées.

Simulation Monte-Carlo

plot of chunk unnamed-chunk-3

Simulation Monte-Carlo

plot of chunk unnamed-chunk-4

Simulation Monte-Carlo

plot of chunk unnamed-chunk-5

Simulation Monte-Carlo

plot of chunk unnamed-chunk-6

Simulation Monte-Carlo

plot of chunk unnamed-chunk-7

Simulation Monte-Carlo

plot of chunk unnamed-chunk-8

Simulation Monte-Carlo

plot of chunk unnamed-chunk-9

Simulation Monte-Carlo

plot of chunk unnamed-chunk-10

Analyse de puissance

plot of chunk unnamed-chunk-11

Calibration de la statistique

Afin de pouvoir comparer plusieurs statistiques ayant subies des contaminations différentes, il est nécessaire de calibrer la statistique :

\[ F=\frac{RSS_{\textrm{petit}}-RSS_{\textrm{grand}}}{RSS_{\textrm{grand}}}\mathbf{\frac{n-p}{p-q}} \]

en \[ F_W=\frac{RSS^{W}_{\textrm{petit}}-RSS^{W}_{\textrm{grand}}}{RSS^{W}_{\textrm{grand}}}\mathbf{\frac{\sum {w_i} -p}{p-q}} \]

Comparaison entre différentes contaminations

plot of chunk unnamed-chunk-12

Comparaison entre différentes contaminations

plot of chunk unnamed-chunk-13

Comparaison entre différentes contaminations

plot of chunk unnamed-chunk-14

Comparaison entre différentes contaminations

plot of chunk unnamed-chunk-15

Conclusion

Produire une version “Kennedy” puis “Huh-Jhun” de la statistique robuste.

Comparer cette statistique au bootstrap robuste.

Comprendre pourquoi utiliser les poids du “grand” modèle ne produit pas un test exact.

Prouver que le test est exact.