Bisakah Anda menghitung kekuatan tes Kolmogorov-Smirnov dalam R?

Apakah mungkin untuk melakukan analisis kekuatan untuk uji Kolmogorov Smirnov 2-sisi pada R?

Saya menguji apakah dua distribusi empiris berbeda menggunakan ks.test (), dan saya ingin menambahkan analisis daya.

Saya belum dapat menemukan analisis daya bawaan untuk tes KS di R. Ada saran?

Sunting : Ini adalah distribusi yang dihasilkan secara acak yang mendekati data saya (dengan ukuran sampel sebenarnya dan perkiraan tingkat peluruhan untuk distribusi eksponensial)

set.seed(100)
x <- rexp(64, rate=0.34)
y <- rexp(54,rate=0.37)

#K-S test: Do x and y come from same distribution?
ks.test(x,y)

Data ini adalah ukuran ukuran tubuh dalam dua kelompok yang berbeda. Saya ingin menunjukkan bahwa kedua kelompok pada dasarnya memiliki distribusi yang sama, tetapi ditanya oleh kolaborator apakah saya memiliki kekuatan untuk mengatakan itu berdasarkan ukuran sampel. Saya secara acak diambil dari distribusi eksponensial di sini, tetapi ini dekat dengan data nyata.

Sejauh ini, saya sudah mengatakan bahwa tidak ada perbedaan yang signifikan dalam distribusi ini berdasarkan pada uji KS dua sisi. Saya juga merencanakan dua distribusi. Bagaimana saya bisa menunjukkan bahwa saya memiliki kekuatan untuk membuat pernyataan seperti itu, mengingat ukuran sampel dan tingkat peluruhan untuk x dan y?

Menemukan kekuatan terhadap alternatif pergeseran skala eksponensial cukup mudah.

Namun, saya tidak tahu bahwa Anda harus menggunakan nilai yang dihitung dari data Anda untuk mengetahui kekuatan apa yang mungkin terjadi. Perhitungan daya post hoc semacam itu cenderung menghasilkan kesimpulan kontra-intuitif (dan mungkin menyesatkan).

Kekuasaan, seperti tingkat signifikansi, adalah fenomena yang Anda hadapi sebelum fakta; Anda akan menggunakan pemahaman apriori (termasuk teori, penalaran atau penelitian sebelumnya) untuk memutuskan serangkaian alternatif yang masuk akal untuk dipertimbangkan, dan ukuran efek yang diinginkan

Anda juga dapat mempertimbangkan berbagai alternatif lain (misalnya Anda dapat menanamkan eksponensial di dalam keluarga gamma untuk mempertimbangkan dampak lebih atau kurang kasus miring).

Pertanyaan biasa yang mungkin bisa dijawab oleh analisis kekuatan adalah:

1) apa kekuatannya, untuk ukuran sampel tertentu, pada beberapa ukuran efek atau serangkaian ukuran efek *?

2) diberi ukuran dan kekuatan sampel, seberapa besar efek yang dapat dideteksi?

3) Diberikan kekuatan yang diinginkan untuk ukuran efek tertentu, ukuran sampel apa yang dibutuhkan?

* (di sini 'ukuran efek' dimaksudkan secara umum, dan mungkin misalnya, rasio rata-rata tertentu, atau perbedaan rata-rata, tidak harus standar).

Jelas Anda sudah memiliki ukuran sampel, jadi Anda tidak dalam kasus (3). Anda mungkin mempertimbangkan kasus (2) atau kasus (1).

Saya menyarankan case (1) (yang juga memberi cara untuk menangani case (2)).

Untuk mengilustrasikan pendekatan pada kasus (1) dan melihat bagaimana hubungannya dengan kasus (2), mari kita pertimbangkan contoh spesifik, dengan:

• alternatif pergeseran skala

• populasi eksponensial

• ukuran sampel dalam dua sampel 64 dan 54

Karena ukuran sampel berbeda, kami harus mempertimbangkan kasus di mana penyebaran relatif di salah satu sampel lebih kecil dan lebih besar dari 1 (jika ukurannya sama, pertimbangan simetri memungkinkan untuk mempertimbangkan hanya satu sisi). Namun, karena mereka cukup dekat dengan ukuran yang sama, efeknya sangat kecil. Bagaimanapun, perbaiki parameter untuk salah satu sampel dan ubah yang lain.

Jadi yang dilakukan adalah:

Sebelumnya:

choose a set of scale multipliers representing different alternatives
select an nsim (say 1000)
set mu1=1

Untuk melakukan perhitungan:

for each possible scale multiplier, kappa 
  repeat nsim times
    generate a sample of size n1 from Exp(mu1) and n2 from Exp(kappa*mu1)
    perform the test
  compute the rejection rate across nsim tests at this kappa

Di R, saya melakukan ini:

alpha = 0.05
n1 = 54
n2 = 64
nsim = 10000
s = c(1.1,1.2,1.5,2,2.5,3) # set up grid for kappa
s = c(1/rev(s),1,s)        #  also below and at 1
rr = array(NA,length(s))   # to hold rejection rates

for(i in seq_along(s)) rr[i]=mean(replicate(nsim,
                                    ks.test(rexp(n1,1),rexp(n2,s[i]))$p.value)<alpha
                                 )

plot(rr~s,log="x",ylim=c(0,1),type="n") #set up plot
points(rr~rev(s),col=3) # plot the reversed case to show the (tiny) asymmetry+noise
points(rr~s,col=1) # plot the "real" case last 
abline(h=alpha,col=8,lty=2) # draw in alpha

yang memberikan kekuatan "kurva" berikut

Sumbu x berada pada skala log, sumbu y adalah tingkat penolakan.

Sulit untuk mengatakan di sini, tetapi titik hitam sedikit lebih tinggi di sebelah kiri daripada di sebelah kanan (yaitu, ada kekuatan fraksional lebih besar ketika sampel yang lebih besar memiliki skala yang lebih kecil).

Menggunakan invers normal cdf sebagai transformasi dari tingkat penolakan, kita dapat membuat hubungan antara tingkat penolakan yang ditransformasikan dan log kappa (kappa ada sdi plot, tetapi sumbu x adalah skala-log) sangat hampir linier (kecuali mendekati 0 ), dan jumlah simulasi cukup tinggi sehingga noise sangat rendah - kita bisa mengabaikannya untuk keperluan saat ini.

Jadi kita bisa menggunakan interpolasi linier. Di bawah ini adalah perkiraan ukuran efek untuk daya 50% dan 80% pada ukuran sampel Anda:

Ukuran efek di sisi lain (grup yang lebih besar memiliki skala lebih kecil) hanya sedikit bergeser dari itu (dapat mengambil ukuran efek yang lebih kecil secara fraksional), tetapi itu membuat sedikit perbedaan, jadi saya tidak akan memaksakan intinya.

Jadi tes akan mengambil perbedaan besar (dari rasio skala 1), tetapi bukan yang kecil.

Sekarang untuk beberapa komentar: Saya tidak berpikir tes hipotesis sangat relevan dengan pertanyaan yang mendasari minat ( apakah mereka cukup mirip? ), Dan akibatnya perhitungan daya ini tidak memberi tahu kami apa pun yang secara langsung relevan dengan pertanyaan itu.

Saya pikir Anda menjawab pertanyaan yang lebih berguna dengan menentukan apa yang menurut Anda "pada dasarnya sama" sebenarnya, secara operasional. Itu - yang diupayakan secara rasional untuk aktivitas statistik - harus mengarah pada analisis data yang bermakna.

Karena Kolmogorov-Smirnov adalah non-parametrik sehingga menurut definisi tidak ada analisis daya yang berlaku. Untuk memiliki semacam estimasi, Anda perlu mengasumsikan model latar belakang (dan dengan demikian mengalihkan dari dunia non-parametrik ...) dan menggunakannya untuk menghitung salah satu dari berikut ini: ukuran sampel, MDE, atau daya (yaitu, Anda perbaiki / pilih dua dan hitung yang ketiga).