Saya menggunakan paket "boot" untuk menghitung nilai p bootstrap 2-sisi yang diperkirakan tetapi hasilnya terlalu jauh dari nilai p menggunakan t.test. Saya tidak tahu apa yang saya lakukan salah dalam kode R. Dapatkah seseorang tolong beri saya petunjuk untuk ini

time = c(14,18,11,13,18,17,21,9,16,17,14,15,
         12,12,14,13,6,18,14,16,10,7,15,10)
group=c(rep(1:2, each=12))
sleep = data.frame(time, group)

require(boot)
diff = function(d1,i){
    d = d1[i,]
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

set.seed(1234)
b3 = boot(data = sleep, statistic = diff, R = 5000, strata=sleep$group)

pvalue = mean(abs(b3$t) > abs(b3$t0))
pvalue 

Nilai p yang di-boot 2 sisi (pvalue) = 0,4804 tetapi nilai p 2 sisi dari t.test adalah 0,04342. Kedua nilai-p ada perbedaan sekitar 11 kali. Bagaimana ini bisa terjadi?

Anda menggunakan bootstrap untuk menghasilkan data di bawah distribusi empiris dari data yang diamati. Ini berguna untuk memberikan interval kepercayaan pada perbedaan antara dua cara:

> quantile(b3$t,c(0.025,0.975))
     2.5%     97.5% 
0.4166667 5.5833333 

Untuk mendapatkan nilai- , Anda perlu membuat permutasi berdasarkan hipotesis nol. Ini bisa dilakukan misalnya seperti ini:$p$

diff2 = function(d1,i){
    d = d1; 
    d$group <- d$group[i];  # randomly re-assign groups
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

> set.seed(1234)
> b4 = boot(data = sleep, statistic = diff2, R = 5000)
> mean(abs(b4$t) > abs(b4$t0))
[1] 0.046

Dalam solusi ini, ukuran grup tidak tetap, Anda secara acak menetapkan ulang grup untuk setiap individu dengan bootstraping dari set grup awal. Tampaknya sah bagi saya, namun solusi yang lebih klasik adalah memperbaiki jumlah individu dari masing-masing kelompok, sehingga Anda hanya mengubah urutan grup daripada bootstraping (ini biasanya dimotivasi oleh desain percobaan, di mana ukuran grup ditetapkan sebelumnya. ):

> R <- 10000; d <- sleep
> b5 <- numeric(R); for(i in 1:R) { 
+    d$group <- sample(d$group, length(d$group)); 
+    b5[i] <- mean(d$time[d$group==1])-mean(d$time[d$group==2]); 
+ }
> mean(abs(b5) > 3)
[1] 0.0372

Jawaban Elvis bergantung pada permutasi tetapi menurut saya tidak menjelaskan apa yang salah dengan pendekatan bootstrap asli. Biarkan saya membahas solusi yang hanya didasarkan pada bootstrap.

Masalah penting dari simulasi asli Anda adalah bahwa bootstrap selalu memberi Anda distribusi statistik uji yang BENAR. Namun, ketika menghitung nilai-p Anda harus membandingkan nilai statistik uji yang diperoleh dengan distribusinya di bawah H0, yaitu tidak dengan distribusi yang sebenarnya!

[Mari kita perjelas. Sebagai contoh, diketahui bahwa statistik uji T dari uji-t klasik memiliki distribusi t "sentral" klasik di bawah H0 dan distribusi noncentral pada umumnya. Namun, semua orang mengetahui fakta bahwa nilai T yang diamati dibandingkan dengan distribusi t "sentral" klasik, yaitu seseorang tidak mencoba untuk mendapatkan distribusi t [noncenral] yang benar untuk membuat perbandingan dengan T.]

Nilai p Anda 0,4804 sangat besar, karena nilai yang diamati "t0" dari statistik uji Mean [1] -Mean [2] terletak sangat dekat dengan pusat sampel bootstrap "t". Itu alami dan biasanya selalu demikian [terlepas dari validitas H0], karena sampel yang di-boot "t" mengemulasi distribusi SEBENARNYA dari Mean [1] -Mean [2]. Tetapi, seperti yang disebutkan di atas [dan juga oleh Elvis], yang Anda perlukan adalah distribusi Mean [1] -Mean [2] DI BAWAH H0. Jelas itu

1) di bawah H0 distribusi Mean [1] -Mean [2] akan berpusat di sekitar 0,

2) bentuknya tidak tergantung pada validitas H0.

Dua poin ini menyiratkan bahwa distribusi Mean [1] -Mean [2] di bawah H0 dapat ditiru oleh sampel bootstrap "t" SHIFTED sehingga berpusat di sekitar 0. Dalam R:

b3.under.H0 <- b3$t - mean(b3$t)

dan nilai-p yang sesuai adalah:

mean(abs(b3.under.H0) > abs(b3$t0))

yang memberi Anda nilai "sangat bagus" dari 0,0232. :-)

Biarkan saya perhatikan bahwa poin "2)" yang disebutkan di atas disebut "translasi kesetaraan" dari statistik uji dan TIDAK harus berlaku secara umum! Yaitu untuk beberapa statistik uji, pemindahan "t" yang di-boot tidak memberi Anda perkiraan yang valid dari distribusi statistik uji dalam HO! Lihatlah diskusi ini dan terutama pada jawaban P. Dalgaard: http://tolstoy.newcastle.edu.au/R/e6/help/09/04/11096.html

Masalah pengujian Anda menghasilkan distribusi simetris yang sempurna dari statistik pengujian, tetapi perlu diingat bahwa ada beberapa masalah dengan mendapatkan nilai-P DUA-SISI dalam kasus distribusi bootstap yang miring pada statistik uji. Sekali lagi, baca tautan di atas.

[Dan akhirnya, saya akan menggunakan tes permutasi "murni" dalam situasi Anda; yaitu bagian kedua dari jawaban Elvis. :-)]

Ada banyak cara untuk menghitung CI bootstrap dan nilai-p. Masalah utama adalah bahwa tidak mungkin bagi bootstrap untuk menghasilkan data di bawah hipotesis nol. Tes permutasi adalah alternatif berbasis resampling yang layak untuk ini. Untuk menggunakan bootstrap yang tepat, Anda harus membuat beberapa asumsi tentang distribusi sampling dari statistik uji.

Sebuah komentar tentang kurangnya invarian pengujian: sangat mungkin untuk menemukan 95% CI tidak termasuk dari null belum ap> 0,05 atau sebaliknya. Untuk mendapatkan persetujuan yang lebih baik, perhitungan sampel bootstrap di bawah nol harus dilakukan dengan daripada . Dengan kata lain jika kerapatan miring tepat di sampel bootstrap, kerapatan harus miring kiri di nol. Sangat tidak mungkin untuk membalikkan tes untuk CI dengan solusi non-analitis (mis. Resampling) seperti ini.β * 0 = β * - $\beta^*_0 = \hat{\beta} - \hat{\beta}^*$$\beta^*_0 = \hat{\beta}^* - \hat{\beta}$

bootstrap normal

Salah satu pendekatan adalah bootstrap normal di mana Anda mengambil mean dan standar deviasi dari distribusi bootstrap, menghitung distribusi sampling di bawah nol dengan menggeser distribusi dan menggunakan persentil normal dari distribusi nol pada titik perkiraan dalam sampel bootstrap asli . Ini adalah pendekatan yang masuk akal ketika distribusi bootstrap normal, inspeksi visual biasanya mencukupi di sini. Hasil menggunakan pendekatan ini biasanya sangat dekat dengan kuat, atau estimasi kesalahan berbasis sandwich yang kuat terhadap heteroskedastisitas dan / atau asumsi varian sampel terbatas. Asumsi statistik uji normal adalah kondisi yang lebih kuat dari asumsi dalam tes bootstrap berikutnya yang akan saya bahas.

bootstrap persentil

Pendekatan lain adalah bootstrap persentil yang menurut saya sebagian besar dari kita pertimbangkan ketika kita berbicara tentang bootstrap. Di sini, distribusi parameter yang bootstrap memperkirakan distribusi empiris sampel di bawah hipotesis alternatif. Distribusi ini mungkin tidak normal. A 95% CI mudah dihitung dengan mengambil kuantil empiris. Tetapi satu asumsi penting adalah bahwa distribusi semacam itu sangat penting . Ini berarti bahwa jika parameter yang mendasarinya berubah, bentuk distribusi hanya bergeser oleh konstanta, dan skala tidak selalu berubah. Ini asumsi yang kuat! Jika ini berlaku, Anda dapat menghasilkan "distribusi statistik di bawah hipotesis nol" (DSNH atau$\mathcal{F}_0^*$) dengan mengurangi distribusi bootstrap dari perkiraan, lalu menghitung berapa persentase DSNH "lebih ekstrem" dari perkiraan Anda dengan menggunakan$2 \times \min (\mathcal{F}_0^*(\hat{\beta}), 1-\mathcal{F}_0^*(\hat{\beta}))$

Bootstrap resmi

Solusi bootstrap termudah untuk menghitung nilai- adalah dengan menggunakan bootstrap mahasiswa. Dengan setiap iterasi bootstrap, hitung statistik dan kesalahan standarnya dan kembalikan statistik siswa. Ini memberikan distribusi siswa bootstrap untuk hipotesis yang dapat digunakan untuk menghitung nilai cis dan p sangat mudah. Ini juga mendasari intuisi di balik bootstrap yang bias dikoreksi. Distribusi t bergeser lebih mudah di bawah nol karena hasil terpencil diturunkan oleh varians tinggi yang sesuai.$p$

Contoh pemrograman

Sebagai contoh, saya akan menggunakan citydata dalam paket bootstrap. Interval kepercayaan bootstrap dihitung dengan kode ini:

ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
city.boot <- boot(city, ratio, R = 999, stype = "w", sim = "ordinary")
boot.ci(city.boot, conf = c(0.90, 0.95),
        type = c("norm", "basic", "perc", "bca"))

dan hasilkan output ini:

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = city.boot, conf = c(0.9, 0.95), type = c("norm", 
    "basic", "perc", "bca"))

Intervals : 
Level      Normal              Basic         
90%   ( 1.111,  1.837 )   ( 1.030,  1.750 )   
95%   ( 1.042,  1.906 )   ( 0.895,  1.790 )  

Level     Percentile            BCa          
90%   ( 1.291,  2.011 )   ( 1.292,  2.023 )   
95%   ( 1.251,  2.146 )   ( 1.255,  2.155 )  
Calculations and Intervals on Original Scale

95% CI untuk bootstrap normal diperoleh dengan menghitung:

with(city.boot, 2*t0 - mean(t) + qnorm(c(0.025, 0.975)) %o% sqrt(var(t)[1,1]))

Nilai p diperoleh demikian:

> with(city.boot, pnorm(abs((2*t0 - mean(t) - 1) / sqrt(var(t)[1,1])), lower.tail=F)*2)
[1] 0.0315

Yang setuju bahwa CI normal 95% tidak termasuk nilai rasio nol 1.

Persentil CI diperoleh (dengan beberapa perbedaan karena metode ikatan):

quantile(city.boot$t, c(0.025, 0.975))

Dan nilai p untuk bootstrap persentil adalah:

cvs <- quantile(city.boot$t0 - city.boot$t + 1, c(0.025, 0.975))
mean(city.boot$t > cvs[1] & city.boot$t < cvs[2])

Memberikan p ap 0,035 yang juga setuju dengan interval kepercayaan dalam hal pengecualian 1 dari nilai. Kita tidak dapat secara umum mengamati bahwa, sementara lebar CI persentil hampir selebar CI normal dan bahwa CI persentil lebih jauh dari nol bahwa CI persentil harus memberikan nilai-p yang lebih rendah. Ini karena bentuk distribusi sampling yang mendasari CI untuk metode persentil tidak normal.