Überraschendes Verhalten der Leistung des Fisher-Exact-Tests (Permutationstests)

9

Ich traf auf ein paradoxes Verhalten von sogenannten "exakten Tests" oder "Permutationstests", deren Prototyp der Fisher-Test ist. Hier ist es.

Stellen Sie sich vor, Sie haben zwei Gruppen von 400 Personen (z. B. 400 Kontrollpersonen gegenüber 400 Fällen) und eine Kovariate mit zwei Modalitäten (z. B. exponiert / nicht exponiert). Es gibt nur 5 exponierte Personen, alle in der zweiten Gruppe. Der Fisher-Test sieht folgendermaßen aus:

> 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
(...)

Aber jetzt gibt es eine gewisse Heterogenität in der zweiten Gruppe (den Fällen), z. B. die Form der Krankheit oder das Rekrutierungszentrum. Es kann in 4 Gruppen von 100 Personen aufgeteilt werden. So etwas wird wahrscheinlich passieren:

> 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
(...)

Jetzt haben wir ...p<0,05

Dies ist nur ein Beispiel. Wir können jedoch die Leistungsfähigkeit der beiden Analysestrategien simulieren, vorausgesetzt, dass bei den ersten 400 Personen die Expositionshäufigkeit 0 und bei den 400 verbleibenden Personen 0,0125 beträgt.

Wir können die Aussagekraft der Analyse mit zwei Gruppen von 400 Personen abschätzen:

> 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

Und mit einer Gruppe von 400 und 4 Gruppen von 100 Personen:

> 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

Es gibt einen ziemlichen Machtunterschied. Die Aufteilung der Fälle in 4 Untergruppen ergibt einen leistungsfähigeren Test, auch wenn es keinen Unterschied in der Verteilung zwischen diesen Untergruppen gibt. Natürlich ist dieser Leistungsgewinn nicht auf eine erhöhte Fehlerrate vom Typ I zurückzuführen.

Ist dieses Phänomen bekannt? Bedeutet das, dass die erste Strategie nicht ausreichend ist? Wäre ein Bootstrap-P-Wert eine bessere Lösung? Alle Ihre Kommentare sind willkommen.

Post Scriptum

Wie von @MartijnWeterings hervorgehoben, liegt ein großer Teil des Grundes für dieses Verhalten (was nicht genau meine Frage ist!) In der Tatsache, dass die wahren Typ-I-Fehler der Schleppanalysestrategien nicht dieselben sind. Dies scheint jedoch nicht alles zu erklären. Ich habe versucht, die ROC-Kurven für zu vergleichen gegen .H 1 : p 0 = 0,05 p 1 = 0,0125H.0::p0=p1=0,005H.1::p0=0,05p1=0,0125

Hier ist mein 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)

Hier ist das Ergebnis:

ROC-Kurven

Wir sehen also, dass ein Vergleich mit demselben wahren Typ-I-Fehler immer noch zu (tatsächlich viel kleineren) Unterschieden führt.

Elvis
quelle
Ich verstehe nicht Die Aufteilung der Fallgruppe kann sinnvoll sein, wenn darin eine gewisse Heterogenität vermutet wird - beispielsweise aus 5 verschiedenen Zentren. Die Aufteilung der "exponierten" Modalität scheint mir nicht sinnvoll zu sein.
Elvis
1
Wenn wir den Unterschied zwischen der ersten und der zweiten Strategie grafisch skizzieren würden. Dann stelle ich mir ein Koordinatensystem mit 5 Achsen (für die Gruppen 400 100 100 100 und 100) mit einem Punkt für die Hypothesenwerte und die Oberfläche vor, der einen Abweichungsabstand darstellt, ab dem die Wahrscheinlichkeit unter einem bestimmten Niveau liegt. Bei der ersten Strategie ist diese Oberfläche ein Zylinder, bei der zweiten Strategie ist diese Oberfläche eine Kugel. Gleiches gilt für die wahren Werte und eine Oberfläche um sie herum für den Fehler. Was wir wollen, ist, dass die Überlappung so klein wie möglich ist.
Sextus Empiricus
1
Ich habe das Ende meiner Frage übernommen und ein wenig mehr Einblick in die Gründe gegeben, warum es einen Unterschied zwischen den beiden Methoden gibt.
Sextus Empiricus
1
Ich glaube, Barnards exakter Test wird verwendet, wenn nur einer der beiden Ränder festgelegt ist. Aber wahrscheinlich werden Sie die gleichen Effekte erhalten.
Sextus Empiricus
1
Eine weitere (interessantere) Bemerkung, die ich machen wollte, ist, dass sich die Leistung tatsächlich verringert, wenn Sie mit p0> p1 testen. Die Leistung steigt also, wenn p1> p0 ist, auf dem gleichen Alpha-Level. Aber die Leistung nimmt ab, wenn p1 <p0 ist (ich bekomme sogar eine Kurve, die unterhalb der Diagonale liegt).
Sextus Empiricus

Antworten:

4

Warum unterscheiden sich die p-Werte?

Es gibt zwei Effekte:

  • χ2

    χ2

    Aus diesem Grund unterscheidet sich der p-Wert um fast den Faktor 2. (nicht genau wegen des nächsten Punktes)

  • Während Sie die 5 0 0 0 0 als ebenso extremen Fall verlieren, erhalten Sie die 1 4 0 0 0 als extremeren Fall als 0 2 1 1 1.

χ2


Codebeispiel:

# 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)

Ausgabe des letzten Bits

> 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

Wie es die Leistung beim Aufteilen von Gruppen beeinflusst

  • Es gibt einige Unterschiede aufgrund der diskreten Schritte in den "verfügbaren" Niveaus der p-Werte und der Konservativität des genauen Tests von Fishers (und diese Unterschiede können ziemlich groß werden).

  • Auch der Fisher-Test passt das (unbekannte) Modell basierend auf den Daten an und verwendet dieses Modell dann zur Berechnung der p-Werte. Das Modell im Beispiel ist, dass es genau 5 exponierte Personen gibt. Wenn Sie die Daten mit einem Binomial für die verschiedenen Gruppen modellieren, erhalten Sie gelegentlich mehr oder weniger als 5 Personen. Wenn Sie den Fischertest darauf anwenden, wird ein Teil des Fehlers angepasst und die Residuen sind im Vergleich zu Tests mit festen Rändern kleiner. Das Ergebnis ist, dass der Test viel zu konservativ und nicht genau ist.

α

H.0H.0H.ein

Es bleibt die Frage, ob dies für alle möglichen Situationen gilt.

3-fache Code-Anpassung Ihrer Leistungsanalyse (und 3 Bilder):

Verwendung von Binomialbeschränkung auf den Fall von 5 exponierten Personen

H.0

Es ist interessant zu sehen, dass der Effekt für den Fall 400-400 (rot) viel stärker ist als für den Fall 400-100-100-100-100 (blau). Daher können wir diese Aufteilung tatsächlich verwenden, um die Leistung zu erhöhen und die Wahrscheinlichkeit zu erhöhen, dass H_0 abgelehnt wird. (Obwohl es uns nicht so wichtig ist, den Fehler vom Typ I wahrscheinlicher zu machen, ist es möglicherweise nicht immer so wichtig, diese Aufteilung vorzunehmen, um die Leistung zu erhöhen.)

unterschiedliche Wahrscheinlichkeiten, H0 abzulehnen

unter Verwendung eines Binomials, das nicht auf 5 exponierte Personen beschränkt ist

Wenn wir ein Binom wie Sie verwenden, ergibt keiner der beiden Fälle 400-400 (rot) oder 400-100-100-100-100 (blau) einen genauen p-Wert. Dies liegt daran, dass der exakte Fisher-Test feste Zeilen- und Spaltensummen voraussetzt, das Binomialmodell jedoch zulässt, dass diese frei sind. Der Fisher-Test passt die Zeilen- und Spaltensummen an und macht den Restterm kleiner als den wahren Fehlerterm.

übermäßig konservativer genauer Fisher-Test

Hat die erhöhte Leistung Kosten?

H.0H.einH.ein

Vergleich von H_0 und 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")

Warum wirkt es sich auf die Leistung aus?

Ich glaube, dass der Schlüssel des Problems in der Differenz der Ergebniswerte liegt, die als "signifikant" ausgewählt wurden. Die Situation besteht aus fünf exponierten Personen, die aus 5 Gruppen mit einer Größe von 400, 100, 100, 100 und 100 gezogen werden. Es können verschiedene Auswahlen getroffen werden, die als "extrem" gelten. Anscheinend steigt die Leistung (selbst wenn der effektive Fehler vom Typ I der gleiche ist), wenn wir uns für die zweite Strategie entscheiden.

Wenn wir den Unterschied zwischen der ersten und der zweiten Strategie grafisch skizzieren würden. Dann stelle ich mir ein Koordinatensystem mit 5 Achsen (für die Gruppen 400 100 100 100 und 100) mit einem Punkt für die Hypothesenwerte und die Oberfläche vor, der einen Abweichungsabstand darstellt, ab dem die Wahrscheinlichkeit unter einem bestimmten Niveau liegt. Bei der ersten Strategie ist diese Oberfläche ein Zylinder, bei der zweiten Strategie ist diese Oberfläche eine Kugel. Gleiches gilt für die wahren Werte und eine Oberfläche um sie herum für den Fehler. Was wir wollen, ist, dass die Überlappung so klein wie möglich ist.

Wir können eine tatsächliche Grafik erstellen, wenn wir ein etwas anderes Problem betrachten (mit geringerer Dimensionalität).

H.0::p=0,5

Beispiel eines Mechanismus

Das Diagramm zeigt, wie die Gruppen von 500 und 500 (anstelle einer einzelnen Gruppe von 1000) verteilt sind.

Der Standardhypothesentest würde (für ein Alpha-Niveau von 95%) beurteilen, ob die Summe von X und Y größer als 531 oder kleiner als 469 ist.

Dies schließt jedoch eine sehr unwahrscheinliche ungleiche Verteilung von X und Y ein.

H.0H.ein

Dies ist jedoch nicht (unbedingt) der Fall, wenn wir die Aufteilung der Gruppen nicht zufällig auswählen und wenn die Gruppen möglicherweise eine Bedeutung haben.

Sextus Empiricus
quelle
Versuchen Sie, meinen Code für die Leistungsschätzung auszuführen, indem Sie nur 0,0125 durch 0,02 ersetzen (entsprechend Ihrem Vorschlag, durchschnittlich 8 exponierte Fälle zu haben): Die 400-gegen-400-Analyse hat eine Leistung von 80%, und die Analyse mit 5 Gruppen hat eine Leistung von 90%.
Elvis
Ich stimme jedoch zu, dass die Statistik in der ersten Situation weniger unterschiedliche Werte annehmen kann und dass dies nicht hilft. Dies reicht jedoch nicht aus, um das Problem zu erklären: Diese Leistungsüberlegenheit kann für alle Fehlerstufen des Typs I beobachtet werden, nicht nur für 0,05. Die Quantile der durch die zweite Strategie erhaltenen p-Werte sind immer niedriger als die durch die erste erhaltenen.
Elvis
Ich glaube, ich stimme dem zu, was Sie sagen. Aber was ist die Schlussfolgerung? Würden Sie empfehlen, die Fallgruppe zufällig in 4 Untergruppen aufzuteilen, um etwas Kraft zu gewinnen? Oder stimmen Sie mir zu, dass dies nicht gerechtfertigt werden kann?
Elvis
Ich denke, das Problem ist nicht, dass der Test mit Fällen, die in 4 Untergruppen aufgeteilt sind, schlechte Eigenschaften haben kann - wir waren uns beide einig, dass sich die Fehlerrate des Typs I gut verhalten sollte. Ich denke, dass das Problem darin besteht, dass der Test mit 400 Kontrollen gegenüber 400 Fällen nicht ausreichend ist. Gibt es eine "saubere" Lösung dafür? Könnte der Bootstrap-P-Wert helfen?
Elvis
(Es tut mir leid, dass meine Frage nicht ganz klar war!)
Elvis