Użycie bootstrapu pod H0 do wykonania testu dla różnicy dwóch środków: zastąpienia w grupach lub w próbce zbiorczej

18

Załóżmy, że mam dane z dwoma niezależnymi grupami:

g1.lengths <- c (112.64, 97.10, 84.18, 106.96, 98.42, 101.66)

g2.lengths <- c (84.44, 82.10, 83.26, 81.02, 81.86, 86.80, 
                     85.84, 97.08, 79.64, 83.32, 91.04, 85.92,
                     73.52, 85.58, 97.70, 89.72, 88.92, 103.72,
                     105.02, 99.48, 89.50, 81.74)

group = rep (c ("g1", "g2"), c (length (g1.lengths), length (g2.lengths)))

lengths = data.frame( lengths = c(g1.lengths, g2.lengths), group)

Oczywiste jest, że wielkość próby na grupę jest tendencyjna, gdy g1 ma 6 obserwacji, a g2 ma 22 . Tradycyjna ANOVA sugeruje, że grupy mają różne środki, gdy wartość krytyczna jest ustawiona na 0,05 (wartość p wynosi 0,0044 ).

summary (aov (lengths~group, data = lengths))  

Biorąc pod uwagę, że moim celem jest porównanie średniej różnicy, takie niezrównoważone i małe próbki danych mogą dawać nieodpowiednie wyniki przy tradycyjnym podejściu. Dlatego chcę wykonać test permutacji i bootstrap.

BADANIE UPRAWNIENIA

Hipoteza zerowa (H0) stwierdza, że ​​średnie grupy są takie same. To założenie w teście permutacji jest uzasadnione poprzez połączenie grup w jedną próbkę. Zapewnia to, że próbki dla dwóch grup zostały pobrane z identycznego rozkładu. Poprzez powtarzanie próbkowania (a ściślej - przetasowanie) z zebranych danych, obserwacje są ponownie przydzielane (tasowane) do próbek w nowy sposób i obliczana jest statystyka testu. Wykonanie tego n razy da rozkład próbkowania statystyk testowych przy założeniu, że H0 ma wartość PRAWDA. Na koniec, pod H0, wartość p jest prawdopodobieństwem, że statystyka testu jest równa lub przekracza wartość obserwowaną.

s.size.g1 <- length (g1.lengths)
s.size.g2 <- length (g2.lengths)

pool <- lengths$lengths
obs.diff.p <- mean (g1.lengths) - mean (g2.lengths)
iterations <- 10000
sampl.dist.p <- NULL

set.seed (5)
for (i in 1 : iterations) {
        resample <- sample (c(1:length (pool)), length(pool))

        g1.perm = pool[resample][1 : s.size.g1]
        g2.perm = pool[resample][(s.size.g1+1) : length(pool)]
        sampl.dist.p[i] = mean (g1.perm) - mean (g2.perm) 
}
p.permute <- (sum (abs (sampl.dist.p) >= abs(obs.diff.p)) + 1)/ (iterations+1)

Podana wartość p testu permutacji wynosi 0,0053 . OK, jeśli zrobiłem to poprawnie, permutacje i parametryczna ANOVA dają prawie identyczne wyniki.

BOOTSTRAP

Przede wszystkim mam świadomość, że bootstrap nie może pomóc, gdy próbki są zbyt małe. Ten post pokazał, że może być jeszcze gorzej i wprowadzać w błąd . Po drugie podkreślono, że test permutacji jest ogólnie lepszy niż bootstrap, gdy głównym celem jest testowanie hipotez. Niemniej jednak ten świetny post rozwiązuje ważne różnice między metodami intensywnie korzystającymi z komputera. Jednak tutaj chciałbym postawić (wierzę) inne pytanie.

Najpierw przedstawię najpopularniejsze podejście do ładowania początkowego (Bootstrap1: ponowne próbkowanie w próbce zbiorczej ):

s.size.g1 <- length (g1.lengths)
s.size.g2 <- length (g2.lengths)

pool <- lengths$lengths
obs.diff.b1 <- mean (g1.lengths) - mean (g2.lengths)
iterations <- 10000
sampl.dist.b1 <- NULL

set.seed (5)
for (i in 1 : iterations) {
        resample <- sample (c(1:length (pool)), length(pool), replace = TRUE) 
        # "replace = TRUE" is the only difference between bootstrap and permutations

        g1.perm = pool[resample][1 : s.size.g1]
        g2.perm = pool[resample][(s.size.g1+1) : length(pool)]
        sampl.dist.b1[i] = mean (g1.perm) - mean (g2.perm) 
}
p.boot1 <- (sum (abs (sampl.dist.b1) >= obs.diff.b1) + 1)/ (iterations+1)

Wartość P wykonywanego w ten sposób ładowania początkowego wynosi 0,005 . Nawet jeśli brzmi to rozsądnie i prawie identycznie jak parametryczna ANOVA i test permutacji, czy właściwe jest uzasadnienie H0 w tym bootstrapie na podstawie tego, że właśnie zebraliśmy próbki, z których pobraliśmy kolejne próbki?

Inne podejście znalazłem w kilku artykułach naukowych. W szczególności widziałem, że badacze modyfikują dane w celu spełnienia H0 przed rozpoczęciem ładowania. Rozglądając się, znalazłem bardzo interesujący post w CV, w którym @ jan.s wyjaśnił niezwykłe wyniki bootstrap w pytaniu, w którym celem było porównanie dwóch średnich. Jednak w tym poście nie opisano, jak wykonać bootstrap, gdy dane są modyfikowane przed bootstrapem. Podejście, w którym dane są modyfikowane przed rozpoczęciem ładowania wygląda następująco:

  1. H0 stwierdza, że ​​średnie dwóch grup są takie same
  2. H0 jest prawdą, jeśli odejmiemy indywidualne obserwacje od średniej próbki zbiorczej

W takim przypadku modyfikacja danych powinna wpływać na średnie grup, a zatem na ich różnicę, ale nie na zmienność w obrębie grup (i między nimi).

  1. Zmodyfikowane dane będą podstawą do dalszego ładowania, z zastrzeżeniami, że próbkowanie jest przeprowadzane w ramach każdej grupy osobno .
  2. Różnicę między wartością początkową g1 i g2 oblicza się i porównuje z zaobserwowaną (niemodyfikowaną) różnicą między grupami.
  3. Proporcja równych lub większej liczby wartości ekstremalnych niż obserwowana podzielona przez liczbę iteracji da wartość p.

Oto kod (Bootstrap2: ponowne próbkowanie w grupach po modyfikacji, że H0 ma wartość PRAWDA ):

s.size.g1 <- length (g1.lengths)
s.size.g2 <- length (g2.lengths)

pool <- lengths$lengths
obs.diff.b2 <- mean (g1.lengths) - mean (g2.lengths)

# make H0 to be true (no difference between means of two groups)
H0 <- pool - mean (pool)

# g1 from H0 
g1.H0 <- H0[1:s.size.g1] 

# g2 from H0
g2.H0 <- H0[(s.size.g1+1):length(pool)]

iterations <- 10000
sampl.dist.b2 <- NULL

set.seed (5)
for (i in 1 : iterations) {
        # Sample with replacement in g1
        g1.boot = sample (g1.H0, replace = T)

        # Sample with replacement in g2
        g2.boot = sample (g2.H0, replace = T)

        # bootstrapped difference
        sampl.dist.b2[i] <- mean (g1.boot) - mean (g2.boot)  
}
p.boot2 <- (sum (abs (sampl.dist.b2) >= obs.diff.b2) + 1)/ (iterations+1)

Tak wykonany bootstrap da wartość p 0,514, która jest ogromnie inna w porównaniu do poprzednich testów. Uważam, że ma to związek z wyjaśnieniem @ jan.s , ale nie mogę zrozumieć, gdzie jest klucz ...

Newbie_R
źródło
1
Ciekawy problem i ładnie zaprezentowany. Bootstrap ma problemy, gdy wielkość próby jest bardzo mała tylko dlatego, że pierwotna próbka nie reprezentuje dobrze populacji. Rozmiar próbki nie musi być zbyt duży, aby bootstrap działał. Twoje próbki wielkości 6 i 22 mogą nie być takie złe. W artykule Efrona (1983) pasek startowy został porównany z formą CV do oszacowania poziomów błędów dla liniowych funkcji dyskryminacyjnych w problemach klasyfikacyjnych z 2 klasami, w których wielkość każdej próbki treningowej jest mniejsza niż 10. Omawiam to w mojej książce Metody ładowania początkowego ( 2007).
Michael R. Chernick,
2
W mojej książce znajduje się także rozdział o niepowodzeniu ładowania początkowego oraz omówienie problemu małej wielkości próbki.
Michael R. Chernick
Pojedynczy wiersz, który trzeba naprawić w bootstrap # 2 podejścia do sprawiają, że praca jest taka: H0 <- pool - mean (pool). Należy go zastąpić H0 <- c(g1.lengths - mean(g1.lengths), g2.lengths - mean(g2.lengths)). Otrzymasz wartość p 0,0023. (To jest to samo, co Zenit wyjaśnił w swojej odpowiedzi.) To wszystko, tylko prosty błąd w kodzie. CC do @MichaelChernick
ameba mówi Przywróć Monikę
czy można te metody obezwładnić? Mam na myśli, czy mogą wykryć JAKĄKOLWIEK różnicę, gdy grupy są bardzo duże: pula> 43k.
Alex Alvarez Pérez

Odpowiedzi:

17

Oto moje zdanie na ten temat, oparte na rozdziale 16 Efrona i Tibshirani's An Introduction to bootstrap (strona 220-224). Krótko mówiąc, twój drugi algorytm ładowania początkowego został nieprawidłowo zaimplementowany, ale ogólna idea jest poprawna.

Podczas przeprowadzania testów ładowania początkowego należy upewnić się, że metoda ponownego próbkowania generuje dane, które odpowiadają hipotezie zerowej. Użyję danych dotyczących snu w R, aby zilustrować ten post. Zauważ, że korzystam ze statystyk z testów uczniowskich, a nie tylko z różnic średnich, co jest zalecane w podręczniku.

Klasyczny test t, który wykorzystuje wynik analityczny do uzyskania informacji o rozkładzie próbkowania statystyki t, daje następujący wynik:

x <- sleep$extra[sleep$group==1] y <- sleep$extra[sleep$group==2]
t.test(x,y)
t = -1.8608, df = 17.776, p-value = 0.07939

n1n2

# pooled sample, assumes equal variance
pooled <- c(x,y)
for (i in 1:10000){
  sample.index <- sample(c(1:length(pooled)),replace=TRUE)
  sample.x <- pooled[sample.index][1:length(x)]
  sample.y <- pooled[sample.index][-c(1:length(y))]
  boot.t[i] <- t.test(sample.x,sample.y)$statistic
}
p.pooled <-  (1 + sum(abs(boot.t) > abs(t.test(x,y)$statistic))) / (10000+1) 
p.pooled
[1] 0.07929207

H0H0H0z¯

x~i=xix¯+z¯
y~i=yiy¯+z¯

x~/y~z¯H0

# sample from H0 separately, no assumption about equal variance
xt <- x - mean(x) + mean(sleep$extra) #
yt <- y - mean(y) + mean(sleep$extra)

boot.t <- c(1:10000)
for (i in 1:10000){
  sample.x <- sample(xt,replace=TRUE)
  sample.y <- sample(yt,replace=TRUE)
  boot.t[i] <- t.test(sample.x,sample.y)$statistic
}
p.h0 <-  (1 + sum(abs(boot.t) > abs(t.test(x,y)$statistic))) / (10000+1)  # 
p.h0
[1] 0.08049195

Tym razem otrzymaliśmy podobne wartości p dla trzech podejść. Mam nadzieję że to pomoże!

Zenit
źródło
1
czy byłbyś uprzejmy i wyjaśnił, dlaczego dodano „1” do następującego: (1 + sum(abs(boot.t) > abs(t.test(x,y)$statistic))) / (10000+1)zamiast czegoś takiego: mean(abs(boot.t) > abs(t.test(x,y)$statistic))Dziękuję za poświęcony czas.
TG_Montana
+1. Czy to podejście bootstrap na zmodyfikowanych danych do próbki z H0 ma określoną nazwę?
ameba mówi Przywróć Monikę
3
H0pvalue=number of times {t>tobs}BB
@amoeba: Nie jestem pewien, czy ta procedura ma formalną lub uzgodnioną nazwę, ale myślę, że można ją opisać jako test ładowania początkowego pod kątem równości środków , a nie dystrybucji . Strona z pełną procedurą brakuje w książkach Google , ale jej motywacja pojawia się na stronie 223. Inny jej opis można znaleźć w tych notatkach, na stronie 13 ( galton.uchicago.edu/~eichler/stat24600/Handouts/bootstrap. pdf ).
Zenit
(+1) Doskonała odpowiedź. Czy mógłbyś wyjaśnić, dlaczego „ten algorytm [ponownie próbkujący same dane bez centrowania] faktycznie sprawdza, czy rozkład xiy są identyczne”? Rozumiem, że ta strategia ponownego próbkowania zapewnia, że rozkłady są takie same, ale dlaczego testuje , że są takie same?
półfinał