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 .$H_0 : p_0 = p_1 = 0.005$$H1 : p_0 = 0.05 \ne p1 = 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:

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

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ść .$\chi^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.$\chi^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”.$\chi^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).$\alpha$

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ę.$H_0$$H_0$$H_a$

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.$H_0$

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)

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.

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.$H_0$$H_a$$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)?$H_0: p=0.5$

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.$H_0$$H_a$

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