Obliczanie wartości p za pomocą bootstrap z R

28

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ć?

Tu.2
źródło
dlaczego b3 $ t0 ma dwa wpisy?
Xi'an
1
to nazwa coli!
Elvis
2
Niepoprawnie obliczasz wartość . Dokumentacja mówi, że jest statystyką obserwowaną, a nie rozkładem zerowym, jak sugerowałaby notacja. Musisz wymyślić oszacowanie odległości próbkowania poniżej zera. Zobacz moją odpowiedź, aby uzyskać więcej informacji. Spróbuj wykonać nieskorygowany test odchylenia. t 0pt0mean(abs(b3$t0) < abs(b3$t-mean(b3$t)))
AdamO,

Odpowiedzi:

31

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
Elvis
źródło
5
Jest to technicznie test permutacji, a nie wartość p bootstrap.
AdamO
@AdamO Zgadzam się, że w tej odpowiedzi przedstawiono test permutacji (i jego nieco zmodyfikowany wariant); dzieje się tak, ponieważ podczas ponownego próbkowania grupy są łączone. Natomiast w teście opartym na bootstrapie wartości dla każdej grupy powinny być próbkowane przy użyciu tylko danych dla tej samej grupy. Oto jedna odpowiedź wyjaśniająca, jak to zrobić: stats.stackexchange.com/a/187630/28666 .
ameba mówi Przywróć Monikę
@amoeba Myślę, że odpowiedź, którą podasz link, jest również testem permutacyjnym, związanym z bootstrapem tylko w takim zakresie, w jakim wymaga ponownego próbkowania. Zgłaszanie jest w porządku, ale do zgłaszania wykorzystywane są dwie metody. Nieparametryczny pasek startowy technicznie nie jest w stanie wygenerować danych przy założeniu zerowej hipotezy. Zobacz moją odpowiedź, aby dowiedzieć się, jak generowane są wartości p z dystrybucji bootstrap.
AdamO
@AdamO Wydaje mi się, że to kwestia terminologii, ale nie rozumiem, jak procedurę opisaną w odpowiedzi na link można nazwać testem „permutacji”, ponieważ nic nie jest tam permutowane: wartości próbek dla każdej grupy są generowane na podstawie danych z tego tylko grupa.
ameba mówi Przywróć Monikę
1
Elvis, myślę, że pierwszy fragment kodu w twojej odpowiedzi to również test permutacyjny. Podczas ponownego próbkowania grupujesz grupy razem! To określa test permutacji.
ameba mówi Przywróć Monikę
25

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

Jan S.
źródło
17

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.

β0=β^β^β0=β^-β^

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

fa02)×min(fa0(β^),1-fa0(β^))

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.

AdamO
źródło
To bardzo interesująca odpowiedź @AdamO, ale czy możesz podać kilka przykładów? Na R możesz użyć funkcji boot.cii użyć argumentu „typ”, aby wybrać studentizowany CI (możesz również wybrać CI BCA). Jak jednak obliczyć wartości p? Czy używasz oszacowania lub statystyki testowej? Miałem podobne pytanie, na które odpowiedź byłaby bardzo wdzięczna.
Kevin Zarca
1
+1 za jasne wyjaśnienie korzyści płynących ze studenckiego paska startowego.
eric_kernfeld
@KevinOunet Podałem dwa przykłady replikacji wartości p z elementów CI w pakiecie rozruchowym. czy to pomaga?
AdamO,
1
Dziękuję @AdamO, to naprawdę pomaga! Czy możesz podać ostatni przykład studenckiego bootstrapu?
Kevin Zarca