Uji apakah variabel mengikuti distribusi yang sama

Jika Anda ingin menguji apakah dua variabel mengikuti distribusi yang sama, apakah akan menjadi tes yang bagus untuk menyortir kedua variabel, dan kemudian memeriksa korelasinya? Jika tinggi (setidaknya 0,9?), Maka variabel yang paling mungkin berasal dari distribusi yang sama.

Dengan distribusi di sini yang saya maksud "normal", "chi-square", "gamma" dll.

Mari kita cari tahu apakah ini ujian yang bagus atau tidak. Ada lebih dari itu hanya dari mengklaim itu buruk atau menunjukkan dalam satu contoh bahwa itu tidak berfungsi dengan baik. Sebagian besar tes bekerja buruk dalam beberapa keadaan, sehingga seringkali kita dihadapkan dengan mengidentifikasi keadaan di mana setiap tes yang diusulkan mungkin menjadi pilihan yang baik.

Deskripsi tes

Seperti halnya tes hipotesis, tes ini terdiri dari (a) hipotesis nol dan alternatif dan (b) statistik uji (koefisien korelasi) yang dimaksudkan untuk membedakan antara hipotesis.

Hipotesis nol adalah bahwa dua variabel berasal dari distribusi yang sama. Tepatnya, mari kita beri nama variabel dan dan menganggap kita memiliki diamati contoh X , disebut x_i = (x_1, x_2, \ ldots, x_ {n_x}) , dan n_y contoh Y , disebut y_i . Hipotesis nol adalah bahwa semua instance X dan Y adalah independen dan terdistribusi secara identik (iid).$X$$Y$$n_x$$X$$x_i = (x_1, x_2, \ldots, x_{n_x})$$n_y$$Y$$y_i$$X$$Y$

Mari kita ambil sebagai hipotesis alternatif bahwa (a) semua instance adalah iid menurut beberapa distribusi mendasarinya dan (b) semua instance adalah iid menurut beberapa distribusi mendasarinya tetapi (c) berbeda dari . (Dengan demikian, kami tidak akan mencari korelasi antara , korelasi antara , korelasi antara dan , atau perbedaan distribusi antara 's atau secara terpisah: yang diasumsikan tidak masuk akal.)$X$$F_X$$Y$$F_Y$$F_X$$F_Y$$x_i$$y_i$$x_i$$y_j$$x$$y$

Statistik uji yang diusulkan mengasumsikan bahwa (sebut nilai umum ini ) dan menghitung koefisien korelasi dari (di mana, seperti biasa, menunjuk terkecil dari data). Sebut ini .$n_x = n_y$$n$[ i ] i th t ( x , y $(x_{[i]}, y_{[i]})$$[i]$$i^\text{th}$$t(x,y)$

Uji permutasi

Dalam situasi ini - tidak peduli apa statistik diusulkan - kita selalu dapat melakukan tes permutasi. Di bawah hipotesis nol, kemungkinan data sama dengan kemungkinan setiap permutasi dari data nilai-nilai. Dengan kata lain, penugasan setengah data ke dan setengah lainnya ke adalah kebetulan acak murni. Ini adalah konsekuensi langsung dan sederhana dari asumsi iid dan hipotesis nol bahwa .( ( x 1 , x 2 , ... , x n ) , ( y 1 , y 2 , ... , y n ) ) 2 n X Y F X = F $t$$\left((x_1, x_2, \ldots, x_n), (y_1, y_2, \ldots, y_n)\right)$$2n$$X$$Y$$F_X=F_Y$

Oleh karena itu, distribusi sampling , tergantung pada pengamatan dan , adalah distribusi semua nilai diperoleh untuk semuapermutasi data. Kami tertarik pada hal ini karena untuk setiap ukuran tes yang dimaksudkan yang dimaksudkan , seperti (sesuai dengan kepercayaan %), kami akan membangun wilayah kritis dua sisi dari distribusi sampel : terdiri dari paling ekstrim % dari nilai yang mungkin darix i y i t ( 2 n ) ! α $t(x,y)$$x_i$$y_i$$t$$(2n)!$$\alpha$95 t 100 α $\alpha = .05$$95$$t$$100\alpha$$t$(di sisi atas, karena korelasi tinggi konsisten dengan distribusi yang sama dan korelasi rendah tidak). Ini adalah bagaimana kita menentukan seberapa besar koefisien korelasi untuk memutuskan data berasal dari distribusi yang berbeda.

Mensimulasikan distribusi sampling nol

Karena(atau, jika Anda suka, , yang menghitung jumlah cara untuk membagi data menjadi dua bagian ukuran ) menjadi besar bahkan untuk kecil , tidak praktis untuk menghitung distribusi sampling secara tepat, jadi kami sampel menggunakan simulasi. (Misalnya, ketika , dan Kira- .) Sekitar seribu sampel sering mencukupi (dan tentu akan untuk eksplorasi yang akan kita lakukan).( 2 $(2n)!$2nnnn=16 ( 2$\binom{2n}{n}/2$$2n$$n$$n$$n=16$(2n)! ≈2,63×10$\binom{2n}{n}/2 = 300\ 540\ 195$$(2n)! \approx 2.63\times 10^{35}$

Ada dua hal yang perlu kita ketahui: pertama, seperti apa distribusi sampling terlihat di bawah hipotesis nol. Kedua, seberapa baik tes ini membedakan antara distribusi yang berbeda?

Ada komplikasi: distribusi pengambilan sampel tergantung pada sifat data. Yang bisa kita lakukan adalah melihat data realistis, dibuat untuk meniru apa pun yang ingin kita pelajari, dan berharap apa yang kita pelajari dari simulasi akan berlaku untuk situasi kita sendiri.

Penerapan

Sebagai ilustrasi, saya telah melakukan pekerjaan ini di R. Itu jatuh secara alami menjadi tiga bagian.

1. Suatu fungsi untuk menghitung statistik uji . Karena saya ingin menjadi sedikit lebih umum, versi saya menangani dataset ukuran yang berbeda ( ) dengan menginterpolasi secara linier di antara nilai-nilai dalam dataset yang lebih besar (diurutkan) untuk membuat kecocokan dengan (diurutkan) dataset yang lebih kecil. Karena ini sudah dilakukan oleh fungsi , saya hanya mengambil hasilnya:$t(x,y)$$n_x \ne n_y$Rqqplot

   test.statistic <- function(x, y) {
     transform <- function(z) -log(1-z^2)/2
     fit <- qqplot(x,y, plot.it=FALSE)
     transform(cor(fit$x, fit$y))
   }

   Sedikit twist - tidak perlu tetapi membantu untuk visualisasi - mengekspresikan kembali koefisien korelasi dengan cara yang akan membuat distribusi statistik nol mendekati simetris. Itu yang transformsedang dilakukan.

2. Simulasi distribusi sampling. Untuk input, fungsi ini menerima jumlah iterasi n.iterbersama dengan dua set data dalam array xdan y. Ini menghasilkan berbagai n.iternilai statistik uji. Cara kerja dalamnya harus transparan, bahkan bagi yang bukan Rpengguna:

   permutation.test <- function(n.iter, x, y) {
     z <- c(x,y)
     n.x <- length(x)
     n.y <- length(y)
     n <- length(z)
     k <- min(n.x, n.y)
     divide <- function() {
       i <- sample.int(n, size=k)
       test.statistic(z[i], z[-i])
     }
     replicate(n.iter, divide())
   }

3. Meskipun hanya itu yang perlu kita lakukan tes, untuk mempelajarinya kita akan ingin mengulangi tes berkali-kali. Jadi, kami melakukan pengujian satu kali dan membungkus kode itu dalam lapisan fungsional ketiga, yang secara umum dinamai di fsini, yang dapat kami panggil berulang kali. Untuk membuatnya cukup umum untuk studi luas, untuk input ia menerima ukuran dataset untuk disimulasikan ( n.xdan n.y), jumlah iterasi untuk setiap tes permutasi ( n.iter), referensi ke fungsi testuntuk menghitung statistik uji (Anda akan melihat untuk sementara mengapa kita mungkin tidak ingin melakukan hard-code ini), dan dua fungsi untuk menghasilkan nilai acak iid, satu untuk ( ) dan satu untuk ( ). Sebuah pilihan$X$dist.x$Y$dist.yplot.it berguna untuk membantu melihat apa yang terjadi.

   f <- function(n.x, n.y, n.iter, test=test.statistic, dist.x=runif, dist.y=runif, 
       plot.it=FALSE) {
     x <- dist.x(n.x)
     y <- dist.y(n.y)
     if(plot.it) qqplot(x,y)
   
     t0 <- test(x,y)
     sim <- permutation.test(n.iter, x, y)
     p <- mean(sim > t0) + mean(sim==t0)/2
     if(plot.it) {
       hist(sim, xlim=c(min(t0, min(sim)), max(t0, max(sim))), 
            main="Permutation distribution")
       abline(v=t0, col="Red", lwd=2)
     }
     return(p)
   }

   Outputnya adalah "nilai-p" yang disimulasikan: proporsi simulasi yang menghasilkan statistik yang terlihat lebih ekstrem daripada yang sebenarnya dihitung untuk data.

Bagian (2) dan (3) sangat umum: Anda dapat melakukan studi seperti ini untuk tes yang berbeda hanya dengan mengganti test.statisticdengan perhitungan lain. Kami melakukannya di bawah.

Hasil pertama

Secara default, kode kami membandingkan data yang diambil dari dua distribusi yang seragam. Saya membiarkannya melakukan itu (untuk , yang merupakan set data yang cukup kecil dan karenanya menyajikan kasus uji yang cukup sulit) dan kemudian mengulanginya untuk perbandingan seragam-normal dan perbandingan seragam-eksponensial. (Distribusi seragam tidak mudah dibedakan dari distribusi normal kecuali jika Anda memiliki sedikit lebih dari nilai, tetapi distribusi eksponensial - memiliki kemiringan tinggi dan ekor kanan panjang - biasanya mudah dibedakan dari distribusi seragam.)$n.x = n.y = 16$$16$

set.seed(17)             # Makes the results reproducible
n.per.rep <- 1000        # Number of iterations to compute each p-value
n.reps <- 1000           # Number of times to call `f`
n.x <- 16; n.y <- 16     # Dataset sizes

par(mfcol=c(2,3))        # Lay results out in three columns
null <- replicate(n.reps, f(n.x, n.y, n.per.rep))
hist(null, breaks=20)
plot(null)

normal <- replicate(n.reps, f(n.x, n.y, n.per.rep, dist.y=rnorm))
hist(normal, breaks=20)
plot(normal)

exponential <- replicate(n.reps, f(n.x, n.y, n.per.rep, dist.y=function(n) rgamma(n, 1)))
hist(exponential, breaks=20)
plot(exponential)

Di sebelah kiri adalah distribusi nol dari nilai-p ketika dan keduanya seragam. Kami berharap histogram mendekati seragam (memberikan perhatian utama pada ujung kiri ekstrim, yang berada dalam kisaran hasil "signifikan") - dan sebenarnya - dan bahwa urutan nilai yang diperoleh selama simulasi, ditunjukkan di bawah ini, terlihat acak - dan itu terjadi. Itu bagus. Ini berarti kita dapat beralih ke langkah selanjutnya untuk mempelajari bagaimana ini berubah ketika dan berasal dari distribusi yang berbeda.$X$$Y$$X$$Y$

Plot tengah menguji varian seragam terhadap normal . Lebih sering daripada tidak, nilai-p lebih rendah dari yang diharapkan. Itu menunjukkan kecenderungan tes ini sebenarnya untuk mendeteksi perbedaan. Tapi itu tidak besar. Sebagai contoh, bilah paling kiri dalam histogram menunjukkan bahwa dari 1000 berjalan (terdiri dari 1000 dataset disimulasikan secara terpisah), nilai-p kurang dari hanya sekitar 110 kali. Jika kita menganggap itu "signifikan," maka tes ini hanya memiliki sekitar % peluang untuk mendeteksi perbedaan antara distribusi seragam dan normal berdasarkan pada$16$$x_i$$16$$y_i$f$0.05$$11$$16$nilai-nilai independen dari masing-masing. Itu daya yang sangat rendah. Tapi mungkin itu tidak bisa dihindari, jadi mari kita lanjutkan.

Plot kanan juga menguji distribusi seragam terhadap eksponensial. Hasil ini aneh. Tes ini cenderung, lebih sering daripada tidak, menyimpulkan bahwa data yang seragam dan data eksponensial terlihat sama. Tampaknya "berpikir" bahwa varian seragam dan eksponensial lebih mirip daripada dua variabel seragam! Apa yang terjadi di sini?

Masalahnya adalah bahwa data dari distribusi eksponensial cenderung memiliki beberapa nilai yang sangat tinggi. Ketika Anda membuat sebar dari mereka terhadap nilai-nilai yang didistribusikan secara seragam, maka akan ada beberapa poin jauh ke kanan atas semua yang lain. Itu sesuai dengan koefisien korelasi yang sangat tinggi. Jadi, setiap kali salah satu distribusi menghasilkan beberapa nilai ekstrem, koefisien korelasi adalah pilihan yang mengerikan untuk mengukur seberapa berbedanya distribusi. Ini mengarah ke masalah yang lebih buruk lagi: ketika ukuran dataset bertambah, kemungkinan mendapatkan beberapa pengamatan ekstrem meningkat. Dengan demikian, kita dapat mengharapkan tes ini berkinerja lebih buruk dan lebih buruk dengan meningkatnya jumlah data. Mengerikan sekali ....

Tes yang lebih baik

Pertanyaan awal telah dijawab dalam negatif. Namun, ada uji kuat yang terkenal untuk membedakan antara distribusi: tes Kolmogorov-Smirnov. Alih-alih koefisien korelasi, itu menghitung deviasi vertikal terbesar dari garis dalam plot QQ mereka. (Ketika data berasal dari distribusi yang sama, plot QQ cenderung mengikuti garis ini. Kalau tidak, itu akan menyimpang di suatu tempat; statistik KS mengambil penyimpangan terbesar seperti itu.)$y=x$

Berikut ini adalah Rimplementasinya:

test.statistic <- function(x, y) {
  ks.test(x,y)$statistic
}

Itu benar: sudah ada di dalam perangkat lunak, jadi kita hanya perlu menyebutnya. Tapi tunggu! Jika Anda membaca manual dengan seksama, Anda akan mengetahui bahwa (a) tes tersebut memberikan nilai-p tetapi (b) bahwa nilai-p (sangat) salah ketika keduanya xdan ymerupakan kumpulan data. Ini dimaksudkan untuk digunakan ketika Anda yakin Anda tahu persis dari mana distribusi data xberasal dan Anda ingin melihat apakah itu benar. Dengan demikian tes tidak mengakomodasi dengan baik ketidakpastian tentang distribusi data yang yberasal.

Tidak masalah! Kerangka uji permutasi masih sama validnya. Dengan membuat perubahan sebelumnya test.statistic, yang harus kita lakukan adalah menjalankan kembali studi sebelumnya, tidak berubah. Inilah hasilnya.

Meskipun distribusi nol tidak seragam (kiri atas), itu cukup seragam di bawah atau lebih, yang mana kita benar-benar peduli dengan nilainya. Pandangan sekilas pada plot di bawahnya (kiri bawah) menunjukkan masalah: statistik KS cenderung mengelompok di sekitar beberapa nilai diskrit. (Masalah ini praktis hilang untuk dataset yang lebih besar.)$p=0.20$

Histogram tengah (seragam vs normal) dan kanan (seragam vs eksponensial) melakukan hal yang benar: pada sebagian besar kasus di mana kedua distribusi berbeda, tes ini menghasilkan nilai-p kecil. Misalnya, ia memiliki peluang % untuk menghasilkan nilai-p kurang dari ketika membandingkan seragam dengan normal berdasarkan pada 16 nilai dari masing-masing. Bandingkan ini dengan piddling % yang dicapai dengan uji koefisien korelasi.$70$$0.05$$11$

Histogram yang tepat tidak begitu baik, tetapi setidaknya itu berada di arah yang benar sekarang! Kami memperkirakan bahwa ia memiliki peluang % untuk mendeteksi perbedaan antara distribusi yang seragam dan eksponensial pada tingkat % dan peluang % untuk membuat deteksi tersebut pada tingkat % (karena dua bilah untuk nilai p kurang dari total lebih dari 500 dari 1000 iterasi).$30$$\alpha=5$$50$$\alpha=10$$0.10$

Kesimpulan

Dengan demikian, masalah dengan uji korelasi bukan karena beberapa kesulitan yang melekat dalam pengaturan ini. Tes korelasi tidak hanya berkinerja sangat buruk, tetapi juga buruk dibandingkan dengan tes yang dikenal luas dan tersedia. (Saya kira itu tidak bisa diterima, artinya akan selalu berkinerja lebih buruk, rata-rata, daripada versi permutasi dari tes KS, menyiratkan tidak ada alasan untuk menggunakannya.)

Tidak, korelasi bukan ujian yang bagus untuk ini.

x <- 1:100 #Uniform
y <- sort(rnorm(100)) #Normal
cor(x,y) #.98

Saya tidak tahu tes yang baik yang membandingkan apakah, misalnya dua distribusi keduanya normal, tetapi mungkin dengan rata-rata berbeda dan sd secara tidak langsung, Anda dapat menguji normalitas masing-masing, secara terpisah, dan jika keduanya tampak normal, tebak mereka berdua.

Jika ada sejumlah variabel yang cukup besar, maka ini mungkin menunjukkan lebih banyak korelasi dengan nilai-nilai yang dipesan. Namun, itu tampaknya bukan metode yang sangat berguna, paling tidak karena memberikan sedikit sarana untuk memperkirakan kepercayaan bahwa mereka mungkin menggunakan model yang sama.

Masalah yang harus Anda alami adalah ketika Anda memiliki model dengan rata-rata dan kemiringan yang serupa, tetapi perbedaan dalam kurtosis, karena jumlah pengukuran yang moderat mungkin cukup cocok untuk terlihat berkorelasi cukup baik.

Tampaknya lebih masuk akal untuk memodelkan kedua variabel terhadap distribusi yang berbeda untuk melihat mana yang paling mungkin untuk masing-masing, dan membandingkan hasilnya.

Mungkin ada beberapa manfaat untuk menormalkan kedua nilai, menyortir dan merencanakan masing-masing - ini akan memungkinkan Anda untuk melihat bagaimana kesesuaian membandingkan - dan Anda dapat merencanakan model yang mungkin untuk keduanya, yang akan terkait dengan apa yang Anda sarankan, tetapi daripada mengharapkan jawaban yang konkret, hanya ide visual tentang kedekatan distribusi.