Zaskakujące zachowanie mocy dokładnego testu Fishera (testy permutacji)

9

Spotkałem paradoksalne zachowanie tak zwanych „testów dokładnych” lub „testów permutacyjnych”, których prototypem jest test Fishera. Oto jest.

Wyobraź sobie, że masz dwie grupy po 400 osobników (np. 400 kontroli i 400 przypadków) oraz zmienną towarzyszącą z dwiema modalnościami (np. Eksponowanymi / nienaświetlonymi). Jest tylko 5 narażonych osób, wszystkie w drugiej grupie. Test Fishera wygląda następująco:

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

Ale teraz w drugiej grupie (przypadkach) występuje pewna niejednorodność, np. Forma choroby lub ośrodek rekrutacyjny. Można go podzielić na 4 grupy po 100 osób. Prawdopodobnie wydarzy się coś takiego:

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

Teraz mamy ...p<0.05

To tylko przykład. Możemy jednak zasymulować moc dwóch strategii analizy, zakładając, że u pierwszych 400 osób częstotliwość ekspozycji wynosi 0, a dla pozostałych 400 osób wynosi 0,0125.

Potęgę analizy możemy oszacować na podstawie dwóch grup po 400 osób:

> 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

I z jedną grupą 400 i 4 grupami po 100 osób:

> 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

Różnica mocy jest dość duża. Podział przypadków na 4 podgrupy daje mocniejszy test, nawet jeśli nie ma różnicy w rozkładzie między tymi podgrupami. Oczywiście tego przyrostu mocy nie można przypisać podwyższonemu poziomowi błędów typu I.

Czy to zjawisko jest dobrze znane? Czy to oznacza, że ​​pierwsza strategia jest słaba? Czy lepsza wartość p byłaby lepszym rozwiązaniem? Wszystkie komentarze są mile widziane.

Post Scriptum

Jak wskazał @MartijnWeterings, duża część przyczyn tego zachowania (co nie jest dokładnie moim pytaniem!) Polega na tym, że prawdziwe błędy typu I strategii analizy holowania nie są takie same. Nie wydaje się to jednak wyjaśniać wszystkiego. Próbowałem porównać krzywe ROC dla vs .H0:p0=p1=0.005H1:p0=0.05p1=0.0125

Oto mój kod.

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)

Oto wynik:

krzywe Roca

Widzimy więc, że porównanie przy tym samym prawdziwym błędzie typu I nadal prowadzi do (rzeczywiście znacznie mniejszych) różnic.

Elvis
źródło
Nie rozumiem. Podział grupy przypadków może mieć sens, gdy podejrzewa się w niej pewną heterogeniczność - powiedzmy, że pochodzą one z 5 różnych centrów. Dzielenie „wyeksponowanej” modalności wydaje mi się nie mieć sensu.
Elvis
1
Gdybyśmy naszkicowali graficznie różnicę między pierwszą i drugą strategią. Następnie wyobrażam sobie układ współrzędnych z 5 osiami (dla grup 400 100 100 100 i 100) z punktem dla wartości hipotez i powierzchni, który przedstawia odległość odchylenia, powyżej której prawdopodobieństwo jest poniżej pewnego poziomu. Przy pierwszej strategii powierzchnia ta jest walcem, przy drugiej strategii powierzchnia ta jest kulą. To samo dotyczy prawdziwych wartości i otaczającej je powierzchni błędu. Chcemy, aby nakładanie się było jak najmniejsze.
Sextus Empiricus
1
Przyjąłem koniec mojego pytania, zapewniając nieco więcej wglądu w uzasadnienie, dlaczego istnieje różnica między tymi dwiema metodami.
Sextus Empiricus
1
Uważam, że dokładny test Barnarda jest stosowany, gdy tylko jeden z dwóch marginesów jest ustalony. Ale prawdopodobnie uzyskasz te same efekty.
Sextus Empiricus
1
Inną (bardziej) interesującą notatką, którą chciałem zrobić, jest fakt, że moc faktycznie zmniejsza się, gdy testujesz z p0> p1. Zatem moc wzrasta, gdy p1> p0, na tym samym poziomie alfa. Ale moc maleje, gdy p1 <p0 (dostaję nawet krzywą poniżej przekątnej).
Sextus Empiricus

Odpowiedzi:

4

Dlaczego wartości p są różne

Działają dwa efekty:

  • Ze względu na dyskrecję wartości wybierasz wektor „najbardziej prawdopodobne” 0 2 1 1 1 wektor. Różniłoby się to jednak od (niemożliwego) 0 1,25 1,25 1,25 1,25, który miałby mniejszą wartość .χ2

    W rezultacie wektor 5 0 0 0 0 nie jest już liczony jako co najmniej ekstremalny przypadek (5 0 0 0 0 ma mniejszy niż 0 2 1 1 1). Tak było wcześniej. Dwustronny testu Fishera na hrabiów stołowych 2x2 oba przypadki 5 ekspozycji będących w pierwszej lub drugiej grupy równie ekstremalne.χ2

    Właśnie dlatego wartość p różni się prawie 2-krotnie (nie dokładnie z powodu następnego punktu)

  • Podczas gdy tracisz 5 0 0 0 0 jako równie skrajny przypadek, zyskujesz 1 4 0 0 0 jako przypadek bardziej ekstremalny niż 0 2 1 1 1.

Różnica polega więc na granicy wartości (lub bezpośrednio obliczonej wartości p stosowanej w implementacji R dokładnego testu Fishera). Jeśli podzielisz grupę 400 na 4 grupy po 100, wówczas różne przypadki będą uważane za bardziej lub mniej „ekstremalne” niż inne. 5 0 0 0 0 jest teraz mniej „ekstremalny” niż 0 2 1 1 1. Ale 1 4 0 0 0 jest bardziej „ekstremalny”.χ2


przykład kodu:

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

wyjście tego ostatniego bitu

> 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

Jak wpływa na moc podczas podziału grup

  • Istnieją pewne różnice z powodu dyskretnych kroków w „dostępnych” poziomach wartości p oraz konserwatywności dokładnego testu Fishera (i różnice te mogą stać się dość duże).

  • również test Fishera pasuje do (nieznanego) modelu opartego na danych, a następnie wykorzystuje ten model do obliczenia wartości p. Model w tym przykładzie jest taki, że jest dokładnie 5 narażonych osób. Jeśli modelujesz dane dwumianem dla różnych grup, od czasu do czasu otrzymasz mniej więcej 5 osób. Kiedy zastosujesz do tego test Fishera, część błędu zostanie dopasowana, a reszty będą mniejsze w porównaniu do testów z ustalonymi marginesami. W rezultacie test jest zbyt konserwatywny, a nie dokładny.

Spodziewałem się, że wpływ na prawdopodobieństwo błędu typu eksperymentu I nie byłby tak wielki, gdybyś losowo podzielił grupy. Jeśli hipoteza zerowa jest prawdziwa, to można spotkać w przybliżeniu procent przypadków znacząca wartość p. W tym przykładzie różnice są duże, jak pokazuje obrazek. Głównym powodem jest to, że przy całkowitych 5 ekspozycjach istnieją tylko trzy poziomy bezwzględnej różnicy (5-0, 4-1, 3-2, 2-3, 1-4, 0-5) i tylko trzy dyskretne p- wartości (w przypadku dwóch grup po 400).α

Najbardziej interesujący jest wykres prawdopodobieństw odrzucenia jeśli jest prawdą i jeśli jest prawdą. W tym przypadku poziom alfa i dyskrecja nie mają tak wielkiego znaczenia (wykreślamy efektywny współczynnik odrzucenia) i wciąż widzimy dużą różnicę.H0H0Ha

Pozostaje pytanie, czy dotyczy to wszystkich możliwych sytuacji.

3-krotna regulacja kodu analizy mocy (i 3 zdjęć):

stosowanie dwumianowe ograniczające do przypadku 5 narażonych osób

Wykresy rzeczywistego prawdopodobieństwa odrzucenia w funkcji wybranej alfa. W dokładnym teście Fishera wiadomo, że wartość p jest dokładnie obliczona, ale występuje tylko kilka poziomów (kroków) tak często, że test może być zbyt konserwatywny w stosunku do wybranego poziomu alfa.H0

Interesujące jest to, że efekt jest znacznie silniejszy w przypadku 400–400 (czerwony) w porównaniu z przypadkiem 400–100–100–100–100 (niebieski). Dlatego możemy rzeczywiście użyć tego podziału, aby zwiększyć moc, aby zwiększyć prawdopodobieństwo odrzucenia H_0. (chociaż nie zależy nam na zwiększeniu prawdopodobieństwa błędu typu I, więc cel podzielenia go w celu zwiększenia mocy nie zawsze musi być tak silny)

różne prawdopodobieństwa odrzucenia H0

stosowanie dwumianu nie ograniczając się do 5 osób narażonych

Jeśli użyjemy dwumianu tak jak ty, wówczas żaden z dwóch przypadków 400-400 (czerwony) lub 400-100-100-100-100 (niebieski) nie podaje dokładnej wartości p. Wynika to z faktu, że dokładny test Fishera zakłada stałe sumy wierszy i kolumn, ale model dwumianowy pozwala na ich swobodę. Test Fishera „dopasuje” sumy wierszy i kolumn, dzięki czemu pozostały okres będzie mniejszy niż rzeczywisty błąd.

zbyt konserwatywny dokładny test Fishera

czy zwiększona moc ma swoją cenę?

Jeśli porównamy prawdopodobieństwo odrzucenia, gdy jest prawdą, a gdy jest prawdą (chcielibyśmy, aby pierwsza wartość była niska, a druga wartość wysoka), wówczas zauważymy, że rzeczywiście moc (odrzucenie, gdy jest prawdą) można zwiększyć bez koszt, który zwiększa błąd typu I.H0HaHa

porównanie H_0 i 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")

Dlaczego wpływa na moc

Uważam, że kluczem do problemu jest różnica wartości wyników, które zostały wybrane jako „znaczące”. Sytuacja polega na tym, że pięć narażonych osobników pochodzi z 5 grup o wielkości 400, 100, 100, 100 i 100. Można dokonać różnych wyborów, które są uważane za „ekstremalne”. najwyraźniej moc wzrasta (nawet gdy efektywny błąd typu I jest taki sam), gdy wybieramy drugą strategię.

Gdybyśmy naszkicowali graficznie różnicę między pierwszą i drugą strategią. Następnie wyobrażam sobie układ współrzędnych z 5 osiami (dla grup 400 100 100 100 i 100) z punktem dla wartości hipotez i powierzchni, który przedstawia odległość odchylenia, powyżej której prawdopodobieństwo jest poniżej pewnego poziomu. Przy pierwszej strategii powierzchnia ta jest walcem, przy drugiej strategii powierzchnia ta jest kulą. To samo dotyczy prawdziwych wartości i otaczającej je powierzchni błędu. Chcemy, aby nakładanie się było jak najmniejsze.

Możemy stworzyć rzeczywistą grafikę, gdy weźmiemy pod uwagę nieco inny problem (o mniejszej wymiarowości).

Wyobraź sobie, że chcemy przetestować proces Bernoulliego , wykonując 1000 eksperymentów. Następnie możemy zastosować tę samą strategię, dzieląc 1000 na grupy na dwie grupy o wielkości 500. Jak to wygląda (niech X i Y będą liczbami w obu grupach)?H0:p=0.5

przykład mechanizmu

Wykres pokazuje, w jaki sposób rozmieszczone są grupy 500 i 500 (zamiast pojedynczej grupy 1000).

Standardowy test hipotez oceniałby (dla 95% poziomu alfa), czy suma X i Y jest większa niż 531 lub mniejsza niż 469.

Ale obejmuje to bardzo mało prawdopodobne nierówne rozmieszczenie X i Y.

Wyobraź sobie przesunięcie rozkładu z na . Wtedy obszary na krawędziach nie mają tak dużego znaczenia, a bardziej okrągła granica byłaby bardziej sensowna.H0Ha

Nie jest to jednak (koniecznie) prawdą, gdy nie wybieramy losowego podziału grup i gdy grupy mogą mieć znaczenie.

Sextus Empiricus
źródło
Spróbuj uruchomić mój kod do oszacowania mocy, po prostu zamieniając 0,0125 na 0,02 (aby dopasować się do twojej sugestii posiadania średnio 8 odsłoniętych przypadków): analiza 400 vs 400 ma moc 80%, a analiza z grupą 5 ma moc 90%.
Elvis
Zgadzam się jednak, że statystyki mogą przyjmować mniej różnych wartości w pierwszej sytuacji i że to nie pomaga. Nie wystarczy to jednak do wyjaśnienia problemu: tę wyższość mocy można zaobserwować dla wszystkich poziomów błędu typu I, nie tylko 0,05. Kwantyle wartości p uzyskane przez drugą strategię są zawsze niższe niż te uzyskane przez pierwszą strategię.
Elvis
Myślę, że zgadzam się z tym, co mówisz. Ale jaki jest wniosek? Czy poleciłbyś losowo podzielić grupę przypadków na 4 podgrupy, aby zyskać trochę mocy? Czy zgadzasz się ze mną, że nie można tego uzasadnić?
Elvis
Myślę, że problemem nie jest to, że test z przypadkami podzielonymi na 4 podgrupy może mieć złe właściwości - obaj zgodziliśmy się co do tego, że poziom błędu typu I powinien zachowywać się dobrze. Myślę, że problem polega na tym, że test z 400 kontrolami w porównaniu do 400 przypadków jest niedociążony. Czy istnieje jakieś „czyste” rozwiązanie tego problemu? Czy wartość p bootstrap może pomóc?
Elvis
(Przykro mi, że moje pytanie nie było w pełni jasne!)
Elvis,