Korzystam z pakietu „boot”, aby obliczyć przybliżoną 2-stronną wartość p ładowania początkowego, ale wynik jest zbyt daleko od wartości p użycia t.test. Nie mogę zrozumieć, co zrobiłem źle w moim kodzie R. Czy ktoś może mi dać na to wskazówkę

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 

2-stronna wartość p ładowania początkowego (wartość) = 0,4804, ale 2-stronna wartość p dla testu t wynosi 0,04342. Obie wartości p różnią się około 11 razy. Jak to się może stać?

Używasz bootstrap do generowania danych w ramach empirycznego rozkładu obserwowanych danych. Może to być przydatne, aby podać przedział ufności dla różnicy między tymi dwoma średnimi:

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

Aby uzyskać wartość , musisz wygenerować permutacje w ramach hipotezy zerowej. Można to zrobić np. W następujący sposób:$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

W tym rozwiązaniu wielkość grup nie jest stała, losowo przypisujesz grupę do każdej osoby, inicjując z początkowego zestawu grup. Wydaje mi się to uzasadnione, jednak bardziej klasycznym rozwiązaniem jest ustalenie liczby osób w każdej grupie, więc po prostu permutuj grupy zamiast bootstrapowania (zwykle jest to motywowane projektem eksperymentu, w którym wielkości grup są wcześniej ustalane ):

> 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

Odpowiedź Elvisa opiera się na permutacjach, ale moim zdaniem nie wyjaśnia, co jest nie tak z oryginalnym podejściem do ładowania początkowego. Pozwól, że omówię rozwiązanie oparte wyłącznie na bootstrapie.

Kluczowym problemem oryginalnej symulacji jest to, że bootstrap zawsze zapewnia PRAWDZIWY rozkład statystyki testowej. Jednak przy obliczaniu wartości p należy porównać uzyskaną wartość statystyki testowej z jej rozkładem WG H0, tzn. Nie z rozkładem prawdziwym!

[Wyjaśnijmy to. Na przykład wiadomo, że statystyka testowa T klasycznego testu t ma klasyczny „centralny” rozkład t pod H0 i ogólnie rozkład niecentralny. Jednak wszyscy są zaznajomieni z faktem, że zaobserwowaną wartość T porównuje się do klasycznego „centralnego” rozkładu t, tj. Nie próbuje się uzyskać prawdziwego [niecenralnego] rozkładu t, aby dokonać porównania z T.]

Twoja wartość p 0,4804 jest tak duża, ponieważ zaobserwowana wartość „t0” statystyki testowej Średnia [1] - Średnia [2] leży bardzo blisko środka próbki „t” ładowanej od początku. Jest to naturalne i zazwyczaj tak jest zawsze [tj. Niezależnie od ważności H0], ponieważ próbka „t” inicjująca ładowanie naśladuje AKTYWNY rozkład średniej [1] - średniej [2]. Ale, jak wspomniano powyżej [a także Elvisa], tak naprawdę potrzebujesz rozkładu średniego [1] - Średniego [2] POD H0. To oczywiste, że

1) pod H0 rozkład średniej [1] - średniej [2] będzie wyśrodkowany wokół 0,

2) jego kształt nie zależy od ważności H0.

Te dwa punkty sugerują, że rozkład średniej [1] - średniej [2] w H0 może być emulowany przez próbkę rozruchową „t” SHIFTED, tak że jest wyśrodkowana wokół 0. W R:

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

a odpowiadająca jej wartość p będzie wynosić:

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

co daje „bardzo ładną” wartość 0,0232. :-)

Pragnę zauważyć, że punkt „2)” wspomniany powyżej nazywa się „ekwiwariantem translacji” statystyki testowej i wcale NIE musi on obowiązywać! Tzn. Dla niektórych statystyk testowych, przesunięcie „t” bootstrapu nie zapewnia prawidłowego oszacowania rozkładu statystyki testowej w HO! Obejrzyj tę dyskusję, a zwłaszcza odpowiedź P. Dalgaarda: http://tolstoy.newcastle.edu.au/R/e6/help/09/04/11096.html

Twój problem z testowaniem daje idealnie symetryczny rozkład statystyki testowej, ale pamiętaj, że istnieją pewne problemy z uzyskaniem DWUSTRONNYCH wartości p w przypadku wypaczonego rozkładu statystyki testowej. Ponownie przeczytaj powyższy link.

[I w końcu użyłbym „czystego” testu permutacji w twojej sytuacji; tj. druga połowa odpowiedzi Elvisa. :-)]

Istnieje wiele sposobów obliczania CI i wartości p bootstrap. Głównym problemem jest to, że bootstrap nie jest w stanie wygenerować danych pod hipotezą zerową. Test permutacji jest realną alternatywą dla tego opartego na ponownym próbkowaniu. Aby użyć prawidłowego paska startowego, należy przyjąć pewne założenia dotyczące rozkładu próbkowania statystyki testowej.

$\beta^*_0 = \hat{\beta} - \hat{\beta}^*$$\beta^*_0 = \hat{\beta}^* - \hat{\beta}$

normalny bootstrap

Jednym z podejść jest normalny bootstrap, w którym bierzesz średnią i standardowe odchylenie rozkładu bootstrap, obliczasz rozkład próbkowania pod wartością zerową, przesuwając rozkład i używając normalnych percentyli z rozkładu zerowego w punkcie oszacowania w oryginalnej próbce bootstrap . Jest to rozsądne podejście, gdy rozkład ładowania jest normalny, zwykle wystarcza tutaj kontrola wzrokowa. Wyniki wykorzystujące to podejście są zwykle bardzo zbliżone do solidnego lub opartego na kanapce oszacowania błędu, który jest odporny na założenia heteroscedastyczności i / lub skończonych wariantów wariancji próbki. Założenie normalnej statystyki testu jest silniejszym warunkiem założeń w następnym teście bootstrapu, który omówię.

percentyl bootstrap

$\mathcal{F}_0^*$$2 \times \min (\mathcal{F}_0^*(\hat{\beta}), 1-\mathcal{F}_0^*(\hat{\beta}))$

Studencki bootstrap

$p$

Przykład programowania

Jako przykład użyję citydanych z pakietu bootstrap. Przedziały ufności ładowania początkowego są obliczane za pomocą tego kodu:

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

i wyprodukuj ten wynik:

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 dla normalnego bootstrap uzyskuje się poprzez obliczenie:

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

W ten sposób uzyskuje się wartość p:

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

Co zgadza się, że 95% normalny CI nie obejmuje wartości zerowej 1.

Uzyskuje się percentyl CI (z pewnymi różnicami wynikającymi z metod wiązania):

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

Wartość p dla bootstrapu percentyla wynosi:

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

Podaje wartość p wynoszącą 0,035, co również zgadza się z przedziałem ufności pod względem wyłączenia 1 z wartości. Zasadniczo nie możemy zaobserwować, że podczas gdy szerokość percentyla CI jest prawie tak szeroka jak normalny CI i że percentyl CI jest dalej od zera, że ​​percentyl CI powinien zapewniać niższe wartości p. Wynika to z faktu, że kształt rozkładu próbkowania leżącego u podstaw CI dla metody percentyla jest nienormalny.