Implementasi tes permutasi mana dalam R untuk digunakan, bukan uji-t (berpasangan dan tidak berpasangan)?

Saya memiliki data dari percobaan yang saya analisis menggunakan uji-t. Variabel dependen adalah skala diskalakan dan data baik tidak berpasangan (yaitu, 2 kelompok) atau berpasangan (yaitu, dalam subyek). Misalnya (dalam mata pelajaran):

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

Namun, data tidak normal sehingga satu reviewer meminta kami untuk menggunakan sesuatu selain dari uji-t. Namun, seperti yang dapat dilihat dengan mudah, data tidak hanya tidak terdistribusi secara normal, tetapi distribusinya tidak sama antara kondisi:

Oleh karena itu, tes nonparametrik yang biasa, Mann-Whitney-U-Test (tidak berpasangan) dan Tes Wilcoxon (berpasangan), tidak dapat digunakan karena memerlukan distribusi yang sama antara kondisi. Oleh karena itu, saya memutuskan bahwa beberapa resampling atau tes permutasi akan menjadi yang terbaik.

Sekarang, saya sedang mencari implementasi R dari yang setara berbasis permutasi dari uji-t, atau saran lain tentang apa yang harus dilakukan dengan data.

Saya tahu bahwa ada beberapa paket-R yang dapat melakukan ini untuk saya (misalnya, koin, perm, exactRankTest, dll.), Tetapi saya tidak tahu mana yang harus dipilih. Jadi, jika seseorang dengan pengalaman menggunakan tes-tes ini bisa memberi saya permulaan, itu akan menjadi ubercool.

PEMBARUAN: Akan ideal jika Anda bisa memberikan contoh bagaimana melaporkan hasil dari tes ini.

Seharusnya tidak terlalu penting karena statistik pengujian akan selalu menjadi perbedaan dalam rata-rata (atau sesuatu yang setara). Perbedaan kecil dapat datang dari penerapan metode Monte-Carlo. Mencoba tiga paket dengan data Anda dengan tes satu sisi untuk dua variabel independen:

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
library(coin)                    # for oneway_test(), pvalue()
pvalue(oneway_test(DV ~ IV, alternative="greater", 
                   distribution=approximate(B=9999)))
[1] 0.00330033

library(perm)                    # for permTS()
permTS(DV ~ IV, alternative="greater", method="exact.mc", 
       control=permControl(nmc=10^4-1))$p.value
[1] 0.003

library(exactRankTests)          # for perm.test()
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.003171822

Untuk memeriksa nilai p yang tepat dengan perhitungan manual dari semua permutasi, saya akan membatasi data ke nilai 9 pertama.

x1 <- x1[1:9]
y1 <- y1[1:9]
DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
pvalue(oneway_test(DV ~ IV, alternative="greater", distribution="exact"))
[1] 0.0945907

permTS(DV ~ IV, alternative="greater", exact=TRUE)$p.value
[1] 0.0945907

# perm.test() gives different result due to rounding of input values
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.1029412

# manual exact permutation test
idx  <- seq(along=DV)                 # indices to permute
idxA <- combn(idx, length(x1))        # all possibilities for different groups

# function to calculate difference in group means given index vector for group A
getDiffM <- function(x) { mean(DV[x]) - mean(DV[!(idx %in% x)]) }
resDM    <- apply(idxA, 2, getDiffM)  # difference in means for all permutations
diffM    <- mean(x1) - mean(y1)       # empirical differencen in group means

# p-value: proportion of group means at least as extreme as observed one
(pVal <- sum(resDM >= diffM) / length(resDM))
[1] 0.0945907

coindan exactRankTestskeduanya dari penulis yang sama, tetapi cointampaknya lebih umum dan luas - juga dalam hal dokumentasi. exactRankTeststidak aktif dikembangkan lagi. Karena itu saya akan memilih coin(juga karena fungsi informatif suka support()), kecuali jika Anda tidak suka berurusan dengan objek S4.

EDIT: untuk dua variabel dependen, sintaksnya adalah

id <- factor(rep(1:length(x1), 2))    # factor for participant
pvalue(oneway_test(DV ~ IV | id, alternative="greater",
                   distribution=approximate(B=9999)))
[1] 0.00810081

Beberapa komentar adalah, saya percaya, secara berurutan.

1) Saya akan mendorong Anda untuk mencoba beberapa tampilan visual dari data Anda, karena mereka dapat menangkap hal-hal yang hilang oleh histogram (seperti grafik), dan saya juga sangat menyarankan Anda memplot pada sumbu berdampingan. Dalam hal ini, saya tidak percaya histogram melakukan pekerjaan yang sangat baik dalam mengkomunikasikan fitur-fitur penting dari data Anda. Sebagai contoh, lihatlah plot-plot box berdampingan:

boxplot(x1, y1, names = c("x1", "y1"))

Atau bahkan stripchart berdampingan:

stripchart(c(x1,y1) ~ rep(1:2, each = 20), method = "jitter", group.names = c("x1","y1"), xlab = "")

Lihatlah pusat, penyebaran, dan bentuk dari ini! Sekitar tiga perempat dari data berada jauh di atas median data . Penyebaran kecil, sedangkan penyebaran sangat besar. Baik dan sangat condong ke kiri, tetapi dengan cara yang berbeda. Sebagai contoh, memiliki lima (!) Nilai berulang nol.y 1 x 1 y 1 x 1 y 1 y $x1$$y1$$x1$$y1$$x1$$y1$$y1$

2) Anda tidak menjelaskan secara detail dari mana data Anda berasal, atau bagaimana mereka diukur, tetapi informasi ini sangat penting ketika tiba saatnya untuk memilih prosedur statistik. Apakah dua sampel Anda di atas independen? Adakah alasan untuk meyakini bahwa distribusi marginal dari kedua sampel harus sama (kecuali untuk perbedaan lokasi, misalnya)? Apa pertimbangan sebelum studi yang membuat Anda mencari bukti perbedaan antara kedua kelompok?

3) Uji-t tidak sesuai untuk data ini karena distribusi marjinal jelas tidak normal, dengan nilai ekstrem di kedua sampel. Jika Anda suka, Anda bisa meminta CLT (karena sampel berukuran sedang) untuk menggunakan uji- (yang akan mirip dengan uji-z untuk sampel besar), tetapi diberi kemiringan (dalam kedua variabel) dari data Anda saya tidak akan menilai banding semacam itu sangat meyakinkan. Tentu, Anda dapat menggunakannya untuk menghitung nilai- , tetapi apa manfaatnya bagi Anda? Jika asumsi tidak terpenuhi maka nilai- hanyalah statistik; itu tidak memberi tahu apa yang Anda (mungkin) ingin tahu: apakah ada bukti bahwa dua sampel berasal dari distribusi yang berbeda.p $z$$p$$p$

4) Tes permutasi juga tidak sesuai untuk data ini. Asumsi tunggal dan sering diabaikan untuk tes permutasi adalah bahwa dua sampel dapat ditukarkan di bawah hipotesis nol. Itu berarti bahwa mereka memiliki distribusi marginal yang identik (di bawah nol). Tetapi Anda berada dalam masalah, karena grafik menunjukkan bahwa distribusi berbeda dalam lokasi dan skala (dan bentuk, juga). Jadi, Anda tidak dapat (secara valid) menguji perbedaan dalam lokasi karena skala berbeda, dan Anda tidak dapat (secara valid) menguji perbedaan dalam skala karena lokasinya berbeda. Ups. Sekali lagi, Anda bisa melakukan tes dan mendapatkan nilai- , tapi lalu bagaimana? Apa yang telah Anda capai?$p$

5) Menurut pendapat saya, data ini adalah contoh sempurna (?) Bahwa gambar yang dipilih dengan baik bernilai 1000 tes hipotesis. Kami tidak perlu statistik untuk membedakan antara pensil dan gudang. Pernyataan yang tepat dalam pandangan saya untuk data ini adalah "Data ini menunjukkan perbedaan yang nyata sehubungan dengan lokasi, skala, dan bentuk." Anda dapat menindaklanjuti dengan statistik deskriptif (kuat) untuk masing-masing untuk mengukur perbedaan, dan menjelaskan apa arti perbedaan dalam konteks studi asli Anda.

6) Peninjau Anda mungkin (dan sayangnya) akan menuntut semacam nilai sebagai prasyarat publikasi. Mendesah! Jika itu saya, mengingat perbedaan sehubungan dengan semua yang saya mungkin akan menggunakan tes Kolmogorov-Smirnov nonparametrik untuk mengeluarkan nilai- yang menunjukkan bahwa distribusi berbeda, dan kemudian melanjutkan dengan statistik deskriptif seperti di atas. Anda perlu menambahkan beberapa noise ke dua sampel untuk menghilangkan ikatan. (Dan tentu saja, ini semua mengasumsikan bahwa sampel Anda independen yang tidak Anda nyatakan secara eksplisit.)$p$$p$

Jawaban ini jauh lebih lama dari yang saya inginkan. Maaf soal itu.

Komentar saya bukan tentang implementasi tes permutasi tetapi tentang masalah yang lebih umum yang diangkat oleh data ini dan diskusi tentang itu, khususnya posting oleh G. Jay Kerns.

Kedua distribusi sebenarnya terlihat sangat mirip dengan saya KECUALI untuk kelompok 0s di Y1, yang jauh berbeda dari pengamatan lain dalam sampel itu (terkecil berikutnya adalah sekitar 50 pada skala 0-100) serta semua yang ada di X1. Saya pertama-tama akan menyelidiki apakah ada yang berbeda dengan pengamatan itu.

Kedua, dengan asumsi 0 itu memang termasuk dalam analisis, mengatakan tes permutasi tidak valid karena distribusi tampaknya berbeda menimbulkan pertanyaan. Jika nol itu benar (distribusi identik), dapatkah Anda (dengan probabilitas yang masuk akal) mendapatkan distribusi yang tampak berbeda seperti keduanya? Menjawab itulah inti dari ujian, bukan? Mungkin dalam kasus ini beberapa orang akan menganggap jawabannya jelas tanpa menjalankan tes, tetapi dengan distribusi yang kecil dan aneh ini, saya tidak berpikir saya akan melakukannya.

Ketika pertanyaan ini muncul lagi, saya dapat menambahkan jawaban lain yang diilhami oleh posting blog baru-baru ini melalui R-Blogger dari Robert Kabacoff, penulis Quick-R dan R in Action menggunakan lmPermpaket.

Namun, metode ini menghasilkan hasil yang sangat kontras (dan sangat tidak stabil) dengan yang dihasilkan oleh coinpaket dalam jawaban @caracakl (nilai p analisis dalam-mata pelajaran adalah 0.008). Analisis mengambil persiapan data dari jawaban @ caracal juga:

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
id <- factor(rep(1:length(x1), 2)) 

library(lmPerm)

summary(aovp( DV ~ IV + Error(id)))

menghasilkan:

> summary(aovp( DV ~ IV + Error(id)))
[1] "Settings:  unique SS "

Error: id
Component 1 :
          Df R Sum Sq R Mean Sq
Residuals 19    15946       839


Error: Within
Component 1 :
          Df R Sum Sq R Mean Sq Iter Pr(Prob)  
IV         1     7924      7924 1004    0.091 .
Residuals 19    21124      1112                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Jika Anda menjalankan ini beberapa kali, nilai-p melompat antara ~ .05 dan ~ .1.

Meskipun ini merupakan jawaban untuk pertanyaan, izinkan saya mengizinkan untuk mengajukan pertanyaan di akhir (saya dapat memindahkan ini ke pertanyaan baru jika diinginkan):
Setiap gagasan mengapa analisis ini sangat tidak stabil dan menghasilkan nilai-p yang sangat berbeda untuk analisis koin? Apakah saya melakukan sesuatu yang salah?

Apakah proporsi skor ini? Jika demikian, Anda tentu tidak boleh menggunakan uji parametrik gaussian, dan sementara Anda bisa melanjutkan dengan pendekatan non-parametrik seperti tes permutasi atau bootstrap sarana, saya sarankan Anda akan mendapatkan lebih banyak kekuatan statistik dengan menggunakan pendekatan parametrik non-gaussian yang cocok. Khususnya, setiap kali Anda dapat menghitung ukuran proporsi dalam unit yang diminati (mis. Peserta dalam percobaan), Anda dapat dan mungkin harus menggunakan model efek campuran yang menentukan pengamatan dengan kesalahan yang terdistribusi secara biner. Lihat Dixon 2004 .

Hanya menambahkan pendekatan lain, ezPermdari ezpaket:

> # preparing the data
> DV <- c(x1, y1)
> IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
> id <- factor(rep(1:length(x1), 2))
> df <- data.frame(id=id,DV=DV,IV=IV)
>
> library(ez)
> ezPerm( data = df, dv = DV, wid = id, within = IV, perms = 1000)
|=========================|100%              Completed after 17 s 
  Effect     p p<.05
1     IV 0.016     *

Hal ini tampaknya konsisten dengan oneway_testsatu coinpaket:

> library(coin)
> pvalue(oneway_test(DV ~ IV | id,  distribution=approximate(B=999999)))
[1] 0.01608002
99 percent confidence interval:
 0.01575782 0.01640682

Namun, perhatikan bahwa ini bukan contoh yang sama yang diberikan oleh @caracal . Dalam contoh, ia termasuk alternative="greater", karena perbedaan p-nilai ~0.008vs ~0.016.

The aovppaket disarankan dalam salah satu jawaban menghasilkan curiga menurunkan p-nilai, dan berjalan curiga cepat bahkan ketika saya mencoba nilai-nilai yang tinggi untuk Iter, Cadan maxIterargumen:

library(lmPerm)
summary(aovp(DV ~ IV + Error(id/IV), data=df,  maxIter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Iter = 1000000000))
summary(aovp(DV ~ IV + Error(id/IV), data=df,  Ca = 0.00000000001))

Yang mengatakan, argumen tampaknya sedikit mengurangi variasi nilai p dari ~.03dan ~.1(saya mendapat kisaran yang lebih besar dari yang dilaporkan oleh @ Henrik), ke 0.03dan 0.07.