Perilaku mengejutkan dari kekuatan uji eksak Fisher (tes permutasi)

9

Saya bertemu dengan perilaku paradoks dari apa yang disebut "tes eksak" atau "tes permutasi", prototipe di antaranya adalah tes Fisher. Ini dia.

Bayangkan Anda memiliki dua kelompok yang terdiri dari 400 orang (mis. 400 kontrol vs 400 kasus), dan kovariat dengan dua modalitas (mis. Terbuka / tidak terpajan). Hanya ada 5 orang yang terpapar, semuanya dalam kelompok kedua. Tes Fisher seperti ini:

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

Tapi sekarang, ada beberapa heterogenitas pada kelompok kedua (kasus), misalnya bentuk penyakit atau pusat perekrutan. Ini dapat dibagi dalam 4 kelompok yang terdiri dari 100 orang. Sesuatu seperti ini kemungkinan akan terjadi:

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

Sekarang, kita memiliki ...p<0.05

Ini hanya sebuah contoh. Tetapi kita dapat mensimulasikan kekuatan dari dua strategi analisis, dengan asumsi bahwa pada 400 individu pertama, frekuensi paparan adalah 0, dan itu adalah 0,0125 pada 400 individu yang tersisa.

Kami dapat memperkirakan kekuatan analisis dengan dua kelompok yang terdiri dari 400 orang:

> 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

Dan dengan satu kelompok yang terdiri dari 400 dan 4 kelompok yang terdiri dari 100 orang:

> 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

Ada perbedaan kekuatan. Membagi kasus dalam 4 subkelompok memberikan tes yang lebih kuat, bahkan jika tidak ada perbedaan distribusi antara subkelompok ini. Tentu saja perolehan daya ini tidak disebabkan oleh peningkatan tingkat kesalahan tipe I.

Apakah fenomena ini terkenal? Apakah itu berarti strategi pertama kurang bertenaga? Apakah nilai p bootstrap menjadi solusi yang lebih baik? Semua komentar Anda dipersilakan.

Skrip Posting

Seperti yang ditunjukkan oleh @ MartijnWeterings, sebagian besar alasan perilaku ini (yang bukan pertanyaan saya!) Terletak pada fakta kesalahan tipe I yang sebenarnya dari strategi analisis derek tidak sama. Namun ini sepertinya tidak menjelaskan semuanya. Saya mencoba membandingkan Kurva ROC untuk vs .H 1 : p 0 = 0,05 p 1 = 0,0125H0:p0=p1=0.005H1:p0=0.05p1=0.0125

Ini kode saya.

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)

Inilah hasilnya:

kurva roc

Jadi kita melihat bahwa perbandingan pada saat yang sama benar kesalahan tipe I masih mengarah ke perbedaan (memang jauh lebih kecil).

Elvis
sumber
Saya tidak mengerti. Memisahkan kelompok kasus dapat masuk akal ketika beberapa heterogenitas diduga ada di dalamnya - katakan, mereka berasal dari 5 pusat yang berbeda. Memisahkan modalitas "terbuka" tampaknya tidak masuk akal bagi saya.
Elvis
1
Jika kita akan membuat sketsa perbedaan antara strategi pertama dan kedua secara grafis. Kemudian saya membayangkan sistem koordinat dengan 5 sumbu (untuk kelompok 400 100 100 100 dan 100) dengan titik untuk nilai-nilai hipotesis dan permukaan yang menggambarkan jarak penyimpangan di luar yang probabilitasnya di bawah tingkat tertentu. Dengan strategi pertama permukaan ini adalah sebuah silinder, dengan strategi kedua permukaan ini adalah sebuah bola. Hal yang sama berlaku untuk nilai-nilai sebenarnya dan permukaan di sekitarnya untuk kesalahan. Yang kami inginkan adalah tumpang tindih sekecil mungkin.
Sextus Empiricus
1
Saya telah mengadopsi akhir dari pertanyaan saya dengan memberikan sedikit lebih banyak wawasan tentang alasan mengapa ada perbedaan antara kedua metode.
Sextus Empiricus
1
Saya percaya tes tepat Barnard digunakan ketika hanya satu dari dua margin yang diperbaiki. Tapi mungkin Anda akan mendapatkan efek yang sama.
Sextus Empiricus
1
satu lagi (lebih) catatan menarik yang ingin saya buat adalah bahwa daya sebenarnya berkurang ketika Anda menguji dengan p0> p1. Jadi daya meningkat ketika p1> p0, pada level alpha yang sama. Tetapi daya berkurang ketika p1 <p0 (saya bahkan mendapatkan kurva di bawah diagonal).
Sextus Empiricus

Jawaban:

4

Mengapa nilai-p berbeda

Ada dua efek yang terjadi:

  • Karena diskritnya nilai Anda memilih vektor 'paling mungkin terjadi' 0 2 1 1 1. Tetapi ini akan berbeda dari (tidak mungkin) 0 1,25 1,25 1,25 1,25, yang akan memiliki nilai lebih kecil .χ2

    Hasilnya adalah bahwa vektor 5 0 0 0 0 tidak dihitung lagi sebagai setidaknya kasus ekstrim (5 0 0 0 0 memiliki dari 0 2 1 1 1). Ini adalah kasus sebelumnya. The dua sisi uji Fisher pada jumlah meja 2x2 kedua kasus dari 5 eksposur berada di pertama atau kelompok kedua sama-sama ekstrem.χ2

    Inilah sebabnya mengapa nilai-p hampir berbeda dengan faktor 2. (bukan karena poin berikutnya)

  • Ketika Anda kehilangan 5 0 0 0 0 sebagai kasus yang sama-sama ekstrim, Anda mendapatkan 1 4 0 0 0 sebagai kasus yang lebih ekstrim daripada 0 2 1 1 1.

Jadi perbedaannya adalah dalam batas nilai (atau nilai-p yang dihitung langsung seperti yang digunakan oleh implementasi R dari uji Fisher yang tepat). Jika Anda membagi kelompok yang terdiri dari 400 menjadi 4 kelompok yang terdiri dari 100, maka kasus yang berbeda akan dianggap lebih atau kurang 'ekstrem' daripada yang lain. 5 0 0 0 0 sekarang kurang 'ekstrem' dari 0 2 1 1 1. Tetapi 1 4 0 0 0 lebih 'ekstrem'.χ2


contoh kode:

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

output dari bit terakhir

> 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

Bagaimana itu mempengaruhi kekuatan ketika memisahkan kelompok

  • Ada beberapa perbedaan karena langkah-langkah yang terpisah dalam tingkat nilai p yang tersedia dan konservatif dari uji eksak Nelayan (dan perbedaan ini mungkin menjadi sangat besar).

  • juga uji Fisher cocok dengan model (tidak diketahui) berdasarkan data dan kemudian menggunakan model ini untuk menghitung nilai-p. Model dalam contoh ini adalah persis ada 5 individu yang terpapar. Jika Anda memodelkan data dengan binomial untuk grup yang berbeda maka Anda akan mendapatkan lebih dari 5 individu. Ketika Anda menerapkan tes fisher untuk ini, maka beberapa kesalahan akan dipasang dan sisanya akan lebih kecil dibandingkan dengan tes dengan marjinal tetap. Hasilnya adalah tes ini terlalu konservatif, tidak tepat.

Saya menduga bahwa efek pada probabilitas kesalahan tipe I eksperimen tidak akan terlalu bagus jika Anda membagi grup secara acak. Jika hipotesis nol itu benar maka Anda akan menemukan dalam persen dari kasus nilai-p yang signifikan. Untuk contoh ini perbedaannya besar seperti yang ditunjukkan gambar. Alasan utama adalah bahwa, dengan total 5 eksposur, hanya ada tiga tingkat perbedaan absolut (5-0, 4-1, 3-2, 2-3, 1-4, 0-5) dan hanya tiga nilai (dalam kasus dua kelompok 400).α

Yang paling menarik adalah plot probabilitas untuk menolak jika benar dan jika benar. Dalam hal ini level alfa dan diskresi tidak terlalu penting (kami merencanakan tingkat penolakan efektif), dan kami masih melihat perbedaan besar.H 0 H aH0H0Ha

Pertanyaannya tetap apakah ini berlaku untuk semua situasi yang mungkin.

3 kali penyesuaian kode analisis daya Anda (dan 3 gambar):

menggunakan pembatasan binomial untuk kasus 5 orang yang terpapar

Plot probabilitas efektif untuk menolak sebagai fungsi dari alpha yang dipilih. Diketahui untuk uji eksak Fisher bahwa nilai-p tepat dihitung tetapi hanya beberapa level (langkah-langkah) yang terjadi sehingga seringkali tes mungkin terlalu konservatif dalam kaitannya dengan tingkat alfa yang dipilih.H0

Sangat menarik untuk melihat bahwa efeknya jauh lebih kuat untuk case 400-400 (merah) dibandingkan case 400-100-100-100-100 (biru). Jadi kita mungkin memang menggunakan pemisahan ini untuk meningkatkan kekuatan, membuatnya lebih cenderung untuk menolak H_0. (walaupun kami tidak terlalu peduli untuk membuat kesalahan tipe I lebih mungkin, jadi titik melakukan pemisahan ini untuk meningkatkan daya mungkin tidak selalu begitu kuat)

probabilitas berbeda untuk menolak H0

menggunakan binomial tidak terbatas pada 5 individu yang terpapar

Jika kami menggunakan binomial seperti yang Anda lakukan maka tak satu pun dari dua kasus 400-400 (merah) atau 400-100-100-100-100 (biru) memberikan nilai-p yang akurat. Ini karena uji eksak Fisher mengasumsikan total baris dan kolom tetap, tetapi model binomial memungkinkan ini menjadi gratis. Uji Fisher akan 'cocok' dengan jumlah baris dan kolom yang membuat istilah residual lebih kecil dari istilah kesalahan sebenarnya.

Tes Fisher yang terlalu konservatif

apakah peningkatan daya datang dengan biaya?

Jika kita membandingkan probabilitas penolakan ketika benar dan ketika benar (kami berharap nilai pertama rendah dan nilai kedua tinggi) maka kita melihat bahwa memang daya (menolak ketika benar) dapat ditingkatkan tanpa biaya bahwa kesalahan tipe I meningkat.H a H aH0HaHa

membandingkan H_0 dan 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")

Mengapa itu mempengaruhi daya

Saya percaya bahwa kunci dari masalah adalah perbedaan nilai hasil yang dipilih untuk menjadi "signifikan". Situasinya adalah lima individu yang terpapar diambil dari 5 kelompok ukuran 400, 100, 100, 100 dan 100. Berbagai pilihan dapat dibuat yang dianggap 'ekstrem'. rupanya kekuatan meningkat (bahkan ketika kesalahan tipe I yang efektif adalah sama) ketika kita memilih strategi kedua.

Jika kita akan membuat sketsa perbedaan antara strategi pertama dan kedua secara grafis. Lalu saya membayangkan sistem koordinat dengan 5 sumbu (untuk kelompok 400 100 100 100 dan 100) dengan titik untuk nilai-nilai hipotesis dan permukaan yang menggambarkan jarak penyimpangan di luar yang probabilitasnya di bawah tingkat tertentu. Dengan strategi pertama permukaan ini adalah sebuah silinder, dengan strategi kedua permukaan ini adalah sebuah bola. Hal yang sama berlaku untuk nilai-nilai sebenarnya dan permukaan di sekitarnya untuk kesalahan. Yang kami inginkan adalah tumpang tindih sekecil mungkin.

Kita dapat membuat grafik aktual ketika kita mempertimbangkan masalah yang sedikit berbeda (dengan dimensi yang lebih rendah).

Bayangkan kita ingin menguji proses Bernoulli dengan melakukan 1000 percobaan. Kemudian kita dapat melakukan strategi yang sama dengan membagi 1000 menjadi kelompok menjadi dua kelompok dengan ukuran 500. Bagaimana ini terlihat (misalkan X dan Y menjadi jumlah di kedua kelompok)?H0:p=0.5

contoh mekanisme

Plot menunjukkan bagaimana kelompok 500 dan 500 (bukan kelompok tunggal 1000) didistribusikan.

Uji hipotesis standar akan menilai (untuk tingkat alfa 95%) apakah jumlah X dan Y lebih besar dari 531 atau lebih kecil dari 469.

Tetapi ini termasuk distribusi X dan Y yang sangat tidak merata.

Bayangkan pergeseran distribusi dari ke . Maka daerah di tepi tidak terlalu penting, dan batas yang lebih melingkar akan lebih masuk akal.H aH0Ha

Namun ini tidak (perlu) benar ketika kita tidak memilih pemisahan kelompok secara acak dan ketika mungkin ada makna untuk kelompok.

Sextus Empiricus
sumber
Cobalah untuk menjalankan kode saya untuk estimasi daya, hanya mengganti 0,0125 dengan 0,02 (untuk mencocokkan saran Anda memiliki rata-rata 8 kasus terbuka): analisis 400 vs 400 memiliki kekuatan 80%, dan analisis dengan 5 grup memiliki kekuatan 90%.
Elvis
Namun saya setuju bahwa statistik dapat mengambil nilai yang kurang berbeda di situasi pertama, dan itu tidak membantu. Namun ini tidak cukup untuk menjelaskan masalah: superioritas daya ini dapat diamati untuk semua tingkat kesalahan tipe I, tidak hanya 0,05. Kuantil dari nilai-p yang diperoleh oleh strategi kedua selalu lebih rendah daripada yang diperoleh oleh yang pertama.
Elvis
Saya pikir saya setuju dengan apa yang Anda katakan. Tapi apa kesimpulannya? Apakah Anda merekomendasikan untuk membagi secara acak kelompok kasus dalam 4 subkelompok, untuk mendapatkan kekuatan? Atau apakah Anda setuju dengan saya bahwa ini tidak dapat dibenarkan?
Elvis
Saya pikir masalahnya bukan bahwa tes dengan kasus terbagi dalam 4 subkelompok mungkin memiliki sifat buruk - kami berdua sepakat pada kenyataan bahwa tingkat kesalahan tipe I harus berperilaku baik. Saya pikir masalahnya adalah bahwa tes dengan 400 kontrol vs 400 kasus kurang bertenaga. Apakah ada solusi "bersih" untuk ini? Bisakah bootstrap nilai p membantu?
Elvis
(Saya minta maaf karena pertanyaan saya tidak sepenuhnya jelas!)
Elvis