Comportement surprenant de la puissance du test exact de Fisher (tests de permutation)

9

J'ai rencontré un comportement paradoxal de soi-disant "tests exacts" ou "tests de permutation", dont le prototype est le test de Fisher. C'est ici.

Imaginez que vous avez deux groupes de 400 individus (par exemple 400 témoins contre 400 cas) et une covariable avec deux modalités (par exemple exposée / non exposée). Il n'y a que 5 individus exposés, tous dans le deuxième groupe. Le test de Fisher ressemble à ceci:

> x <- matrix( c(400, 395, 0, 5) , ncol = 2)
> x
     [,1] [,2]
[1,]  400    0
[2,]  395    5
> fisher.test(x)

    Fisher's Exact Test for Count Data

data:  x
p-value = 0.06172
(...)

Mais maintenant, il y a une certaine hétérogénéité dans le deuxième groupe (les cas), par exemple la forme de la maladie ou le centre de recrutement. Il peut être divisé en 4 groupes de 100 individus. Quelque chose comme ça est susceptible de se produire:

> x <- matrix( c(400, 99, 99 , 99, 98, 0, 1, 1, 1, 2) , ncol = 2)
> x
     [,1] [,2]
[1,]  400    0
[2,]   99    1
[3,]   99    1
[4,]   99    1
[5,]   98    2
> fisher.test(x)

    Fisher's Exact Test for Count Data

data:  x 
p-value = 0.03319
alternative hypothesis: two.sided
(...)

Maintenant, nous avons ...p<0,05

Ce n'est qu'un exemple. Mais nous pouvons simuler la puissance des deux stratégies d'analyse, en supposant que chez les 400 premiers individus, la fréquence d'exposition est de 0, et qu'elle est de 0,0125 chez les 400 individus restants.

On peut estimer la puissance de l'analyse avec deux groupes de 400 individus:

> p1 <- replicate(1000, { n <- rbinom(1, 400, 0.0125); 
                          x <- matrix( c(400, 400 - n, 0, n), ncol = 2); 
                          fisher.test(x)$p.value} )
> mean(p1 < 0.05)
[1] 0.372

Et avec un groupe de 400 et 4 groupes de 100 personnes:

> p2 <- replicate(1000, { n <- rbinom(4, 100, 0.0125); 
                          x <- matrix( c(400, 100 - n, 0, n), ncol = 2);
                          fisher.test(x)$p.value} )
> mean(p2 < 0.05)
[1] 0.629

Il y a une grande différence de puissance. La division des cas en 4 sous-groupes donne un test plus puissant, même s'il n'y a pas de différence de distribution entre ces sous-groupes. Bien entendu, ce gain de puissance n'est pas attribuable à une augmentation du taux d'erreur de type I.

Ce phénomène est-il bien connu? Est-ce à dire que la première stratégie est sous-alimentée? Une valeur p amorcée serait-elle une meilleure solution? Tous vos commentaires sont les bienvenus.

Post Scriptum

Comme l'a souligné @MartijnWeterings, une grande partie de la raison de ce comportement (ce qui n'est pas exactement ma question!) Réside dans le fait que les véritables erreurs de type I des stratégies d'analyse de remorquage ne sont pas les mêmes. Mais cela ne semble pas tout expliquer. J'ai essayé de comparer les courbes ROC pour vs .H 1 : p 0 = 0,05 p 1 = 0,0125H0:p0=p1=0,005H1:p0=0,05p1=0,0125

Voici mon code.

B <- 1e5
p0 <- 0.005
p1 <- 0.0125

# simulation under H0 with p = p0 = 0.005 in all groups
# a = 2 groups 400:400, b = 5 groupe 400:100:100:100:100

p.H0.a <- replicate(B, { n <- rbinom( 2, c(400,400), p0);
                           x <- matrix( c( c(400,400) -n, n ), ncol = 2);
                          fisher.test(x)$p.value} )

p.H0.b <- replicate(B, { n <- rbinom( 5, c(400,rep(100,4)), p0);
                           x <- matrix( c( c(400,rep(100,4)) -n, n ), ncol = 2);
                          fisher.test(x)$p.value} )

# simulation under H1 with p0 = 0.005 (controls) and p1 = 0.0125 (cases)

p.H1.a <- replicate(B, { n <- rbinom( 2, c(400,400), c(p0,p1) );
                           x <- matrix( c( c(400,400) -n, n ), ncol = 2);
                          fisher.test(x)$p.value} )

p.H1.b <- replicate(B, { n <- rbinom( 5, c(400,rep(100,4)), c(p0,rep(p1,4)) );
                           x <- matrix( c( c(400,rep(100,4)) -n, n ), ncol = 2);
                          fisher.test(x)$p.value} )

# roc curve 

ROC <- function(p.H0, p.H1) {
  p.threshold <- seq(0, 1.001, length=501)
  alpha <- sapply(p.threshold, function(th) mean(p.H0 <= th) )
  power <- sapply(p.threshold, function(th) mean(p.H1 <= th) )
  list(x = alpha, y = power)
}

par(mfrow=c(1,2))
plot( ROC(p.H0.a, p.H1.a) , type="b", xlab = "alpha", ylab = "1-beta" , xlim=c(0,1), ylim=c(0,1), asp = 1)
lines( ROC(p.H0.b, p.H1.b) , col="red", type="b" )
abline(0,1)

plot( ROC(p.H0.a, p.H1.a) , type="b", xlab = "alpha", ylab = "1-beta" , xlim=c(0,.1) )
lines( ROC(p.H0.b, p.H1.b) , col="red", type="b" )
abline(0,1)

Voici le résultat:

courbes roc

Nous voyons donc qu'une comparaison avec la même erreur de type I véritable conduit toujours à des différences (bien plus faibles).

Elvis
la source
Je ne comprends pas. La division du groupe de cas peut avoir du sens quand une hétérogénéité est suspectée à l'intérieur - disons, ils proviennent de 5 centres différents. Le fractionnement de la modalité «exposée» ne me semble pas logique.
Elvis
1
Si nous esquissions graphiquement la différence entre la première et la deuxième stratégie. Ensuite, j'imagine un système de coordonnées à 5 axes (pour les groupes de 400100100100 et 100) avec un point pour les valeurs d'hypothèse et la surface qui représente une distance de déviation au-delà de laquelle la probabilité est inférieure à un certain niveau. Avec la première stratégie, cette surface est un cylindre, avec la deuxième stratégie, cette surface est une sphère. La même chose est vraie pour les vraies valeurs et une surface autour d'elle pour l'erreur. Ce que nous voulons, c'est que le chevauchement soit le plus faible possible.
Sextus Empiricus
1
J'ai adopté la fin de ma question en fournissant un peu plus de détails sur le raisonnement pour lequel il y a une différence entre les deux méthodes.
Sextus Empiricus
1
Je crois que le test exact de Barnard est utilisé lorsqu'une seule des deux marges est fixée. Mais vous obtiendrez probablement les mêmes effets.
Sextus Empiricus
1
une autre note (plus) intéressante que je voulais faire est que la puissance diminue réellement lorsque vous testez avec p0> p1. La puissance augmente donc lorsque p1> p0, au même niveau alpha. Mais la puissance diminue lorsque p1 <p0 (j'obtiens même une courbe en dessous de la diagonale).
Sextus Empiricus

Réponses:

4

Pourquoi les valeurs de p sont-elles différentes

Il y a deux effets:

  • χ2

    χ2

    C'est pourquoi la valeur de p diffère de près d'un facteur 2. (pas exactement à cause du point suivant)

  • Pendant que vous perdez le 5 0 0 0 0 comme un cas tout aussi extrême, vous gagnez le 1 4 0 0 0 comme un cas plus extrême que le 0 2 1 1 1.

χ2


exemple de code:

# probability of distribution a and b exposures among 2 groups of 400
draw2 <- function(a,b) {
  choose(400,a)*choose(400,b)/choose(800,5)
}

# probability of distribution a, b, c, d and e exposures among 5 groups of resp 400, 100, 100, 100, 100
draw5 <- function(a,b,c,d,e) {
choose(400,a)*choose(100,b)*choose(100,c)*choose(100,d)*choose(100,e)/choose(800,5)
}

# looping all possible distributions of 5 exposers among 5 groups
# summing the probability when it's p-value is smaller or equal to the observed value 0 2 1 1 1
sumx <- 0
for (f in c(0:5)) {
  for(g in c(0:(5-f))) {
    for(h in c(0:(5-f-g))) {
      for(i in c(0:(5-f-g-h))) {
        j = 5-f-g-h-i
        if (draw5(f, g, h, i, j) <= draw5(0, 2, 1, 1, 1)) {
          sumx <- sumx + draw5(f, g, h, i, j)
        }
      }
    }
  } 
}
sumx  #output is 0.3318617

# the split up case (5 groups, 400 100 100 100 100) can be calculated manually
# as a sum of probabilities for cases 0 5 and 1 4 0 0 0 (0 5 includes all cases 1 a b c d with the sum of the latter four equal to 5)
fisher.test(matrix( c(400, 98, 99 , 99, 99, 0, 2, 1, 1, 1) , ncol = 2))[1]
draw2(0,5) + 4*draw(1,4,0,0,0)

# the original case of 2 groups (400 400) can be calculated manually
# as a sum of probabilities for the cases 0 5 and 5 0 
fisher.test(matrix( c(400, 395, 0, 5) , ncol = 2))[1]
draw2(0,5) + draw2(5,0)

sortie de ce dernier bit

> fisher.test(matrix( c(400, 98, 99 , 99, 99, 0, 2, 1, 1, 1) , ncol = 2))[1]
$p.value
[1] 0.03318617

> draw2(0,5) + 4*draw(1,4,0,0,0)
[1] 0.03318617

> fisher.test(matrix( c(400, 395, 0, 5) , ncol = 2))[1]
$p.value
[1] 0.06171924

> draw2(0,5) + draw2(5,0)
[1] 0.06171924

Comment cela affecte la puissance lors de la division des groupes

  • Il existe certaines différences en raison des étapes discrètes des niveaux «disponibles» des valeurs de p et de la prudence du test exact de Fishers (et ces différences peuvent devenir assez importantes).

  • le test de Fisher correspond également au modèle (inconnu) basé sur les données, puis utilise ce modèle pour calculer les valeurs de p. Le modèle de l'exemple est qu'il y a exactement 5 individus exposés. Si vous modélisez les données avec un binôme pour les différents groupes, vous obtiendrez parfois plus ou moins de 5 individus. Lorsque vous appliquez le test de Fisher à cela, une partie de l'erreur sera ajustée et les résidus seront plus petits par rapport aux tests avec des marginaux fixes. Le résultat est que le test est beaucoup trop conservateur, pas exact.

α

H0H0Hune

Reste à savoir si cela vaut pour toutes les situations possibles.

Ajustement du code 3 fois de votre analyse de puissance (et 3 images):

utilisation d'un binôme restreint au cas de 5 individus exposés

H0

Il est intéressant de voir que l'effet est beaucoup plus fort pour le boîtier 400-400 (rouge) que pour le boîtier 400-100-100-100-100 (bleu). Ainsi, nous pouvons en effet utiliser cette séparation pour augmenter la puissance, la rendre plus susceptible de rejeter le H_0. (bien que nous ne nous soucions pas tellement de rendre l'erreur de type I plus probable, le point de faire cette division pour augmenter la puissance n'est pas toujours si fort)

différentes probabilités de rejeter H0

utilisant un binôme ne se limitant pas à 5 individus exposés

Si nous utilisons un binôme comme vous l'avez fait, aucun des deux cas 400-400 (rouge) ou 400-100-100-100-100 (bleu) ne donne une valeur de p précise. En effet, le test exact de Fisher suppose des totaux de ligne et de colonne fixes, mais le modèle binomial permet de les libérer. Le test de Fisher «ajustera» les totaux des lignes et des colonnes, ce qui rend le terme résiduel plus petit que le vrai terme d'erreur.

test exact de Fisher trop conservateur

l'augmentation de puissance a-t-elle un coût?

H0HuneHune

comparer H_0 et H_a

# using binomial distribution for 400, 100, 100, 100, 100
# x uses separate cases
# y uses the sum of the 100 groups
p <- replicate(4000, { n <- rbinom(4, 100, 0.006125); m <- rbinom(1, 400, 0.006125); 
x <- matrix( c(400 - m, 100 - n, m, n), ncol = 2);
y <- matrix( c(400 - m, 400 - sum(n), m, sum(n)), ncol = 2);
c(sum(n,m),fisher.test(x)$p.value,fisher.test(y)$p.value)} )

# calculate hypothesis test using only tables with sum of 5 for the 1st row
ps <- c(1:1000)/1000
m1 <- sapply(ps,FUN = function(x) mean(p[2,p[1,]==5] < x))
m2 <- sapply(ps,FUN = function(x) mean(p[3,p[1,]==5] < x))

plot(ps,ps,type="l",
     xlab = "chosen alpha level",
     ylab = "p rejection")
lines(ps,m1,col=4)
lines(ps,m2,col=2)

title("due to concervative test p-value will be smaller\n leading to differences")

# using all samples also when the sum exposed individuals is not 5
ps <- c(1:1000)/1000
m1 <- sapply(ps,FUN = function(x) mean(p[2,] < x))
m2 <- sapply(ps,FUN = function(x) mean(p[3,] < x))

plot(ps,ps,type="l", 
     xlab = "chosen alpha level",
     ylab = "p rejection")
lines(ps,m1,col=4)
lines(ps,m2,col=2)

title("overly conservative, low effective p-values \n fitting marginals makes residuals smaller than real error")


#   
# Third graph comparing H_0 and H_a
#
# using binomial distribution for 400, 100, 100, 100, 100
# x uses separate cases
# y uses the sum of the 100 groups
offset <- 0.5
p <- replicate(10000, { n <- rbinom(4, 100, offset*0.0125); m <- rbinom(1, 400, (1-offset)*0.0125); 
x <- matrix( c(400 - m, 100 - n, m, n), ncol = 2);
y <- matrix( c(400 - m, 400 - sum(n), m, sum(n)), ncol = 2);
c(sum(n,m),fisher.test(x)$p.value,fisher.test(y)$p.value)} )

# calculate hypothesis test using only tables with sum of 5 for the 1st row
ps <- c(1:10000)/10000
m1 <- sapply(ps,FUN = function(x) mean(p[2,p[1,]==5] < x))
m2 <- sapply(ps,FUN = function(x) mean(p[3,p[1,]==5] < x))

offset <- 0.6
p <- replicate(10000, { n <- rbinom(4, 100, offset*0.0125); m <- rbinom(1, 400, (1-offset)*0.0125); 
x <- matrix( c(400 - m, 100 - n, m, n), ncol = 2);
y <- matrix( c(400 - m, 400 - sum(n), m, sum(n)), ncol = 2);
c(sum(n,m),fisher.test(x)$p.value,fisher.test(y)$p.value)} )

# calculate hypothesis test using only tables with sum of 5 for the 1st row
ps <- c(1:10000)/10000
m1a <- sapply(ps,FUN = function(x) mean(p[2,p[1,]==5] < x))
m2a <- sapply(ps,FUN = function(x) mean(p[3,p[1,]==5] < x))

plot(ps,ps,type="l",
     xlab = "p rejecting if H_0 true",
     ylab = "p rejecting if H_a true",log="xy")
points(m1,m1a,col=4)
points(m2,m2a,col=2)

legend(0.01,0.001,c("400-400","400-100-100-100-100"),pch=c(1,1),col=c(2,4))

title("comparing H_0:p=0.5 \n with H_a:p=0.6")

Pourquoi cela affecte-t-il la puissance

Je pense que la clé du problème réside dans la différence des valeurs de résultat qui sont choisies comme "significatives". La situation est de cinq individus exposés tirés de 5 groupes de taille 400, 100, 100, 100 et 100. Différentes sélections peuvent être considérées comme «extrêmes». apparemment, la puissance augmente (même lorsque l'erreur efficace de type I est la même) lorsque nous optons pour la deuxième stratégie.

Si nous esquissions graphiquement la différence entre la première et la deuxième stratégie. Ensuite, j'imagine un système de coordonnées à 5 axes (pour les groupes de 400100100100 et 100) avec un point pour les valeurs d'hypothèse et la surface qui représente une distance de déviation au-delà de laquelle la probabilité est inférieure à un certain niveau. Avec la première stratégie, cette surface est un cylindre, avec la deuxième stratégie, cette surface est une sphère. La même chose est vraie pour les vraies valeurs et une surface autour d'elle pour l'erreur. Ce que nous voulons, c'est que le chevauchement soit le plus faible possible.

Nous pouvons faire un graphique réel lorsque nous considérons un problème légèrement différent (avec une dimensionnalité inférieure).

H0:p=0,5

exemple de mécanisme

Le graphique montre comment les groupes de 500 et 500 (au lieu d'un seul groupe de 1000) sont distribués.

Le test d'hypothèse standard évaluerait (pour un niveau alpha de 95%) si la somme de X et Y est supérieure à 531 ou inférieure à 469.

Mais cela inclut une distribution inégale très improbable de X et Y.

H0Hune

Ce n'est cependant pas (nécésarilly) vrai lorsque nous ne sélectionnons pas la division des groupes au hasard et lorsqu'il peut y avoir un sens aux groupes.

Sextus Empiricus
la source
Essayez d'exécuter mon code pour l'estimation de la puissance, en remplaçant simplement 0,0125 par 0,02 (pour correspondre à votre suggestion d'avoir une moyenne de 8 cas exposés): l'analyse 400 vs 400 a une puissance de 80%, et l'analyse avec le groupe 5 a une puissance de 90%.
Elvis
Cependant, je suis d'accord que la statistique peut prendre des valeurs moins différentes dans la première situation, et qu'elle n'aide pas. Cependant, cela ne suffit pas à expliquer le problème: cette supériorité de puissance peut être observée pour tous les niveaux d'erreur de type I, pas seulement 0,05. Les quantiles des valeurs de p obtenues par la deuxième stratégie sont toujours inférieurs à ceux obtenus par la première.
Elvis
Je pense que je suis d'accord avec ce que vous dites. Mais quelle conclusion? Recommanderiez-vous de diviser au hasard le groupe de cas en 4 sous-groupes, pour gagner en puissance? Ou êtes-vous d'accord avec moi que cela ne peut pas être justifié?
Elvis
Je pense que le problème n'est pas que le test avec des cas divisés en 4 sous-groupes peut avoir de mauvaises propriétés - nous nous sommes tous deux entendus sur le fait que son taux d'erreur de type I devrait bien se comporter. Je pense que le problème est que le test avec 400 contrôles contre 400 cas est sous-alimenté. Y a-t-il une solution «propre» à cela? Bootstrap p-value pourrait-il aider?
Elvis
(Je suis désolé que ma question n'ait pas été complètement claire!)
Elvis