Dlaczego testy hipotez na ponownie próbkowanych zestawach danych zbyt często odrzucają wartość zerową?

10

tl; dr: Zaczynając od zestawu danych wygenerowanego pod wartością zerową, ponownie próbkowałem przypadki z zamianą i przeprowadzałem test hipotez na każdym zestawie danych ponownie próbkowanym. Te testy hipotez odrzucają zero w ponad 5% przypadków.

W poniższej, bardzo prostej symulacji, zestawy danych za pomocą $X \sim N(0,1) \amalg Y \sim N(0,1)$ i dopasowuję do nich prosty model OLS. Następnie dla każdego zestawu danych generuję 1000 nowych zestawów danych poprzez ponowne próbkowanie wierszy oryginalnego zestawu danych z zastąpieniem (algorytm opisany szczegółowo w klasycznym tekście Davisona i Hinkleya jako odpowiedni dla regresji liniowej). Do każdego z nich pasuję do tego samego modelu OLS. Ostatecznie około 16% testów hipotez w próbkach bootstrap odrzuca wartość zerową , podczas gdy powinniśmy otrzymać 5% (tak jak robimy to w oryginalnych zestawach danych).

Podejrzewałem, że ma to coś wspólnego z powtarzającymi się obserwacjami powodującymi zawyżone skojarzenia, więc dla porównania wypróbowałem dwa inne podejścia w poniższym kodzie (skomentowałem). W metodzie 2 naprawiam $X$ , a następnie zamieniam $Y$ resztki próbek z modelu OLS w oryginalnym zestawie danych. W metodzie 3 rysuję losową podpróbkę bez zamiany. Obie te alternatywy działają, tj. Ich testy hipotez odrzucają zerowy 5% czasu.

Moje pytanie: czy mam rację, że winowajcami są powtarzające się obserwacje? Jeśli tak, biorąc pod uwagę, że jest to standardowe podejście do ładowania, to gdzie dokładnie naruszamy standardową teorię ładowania?

Aktualizacja nr 1: Więcej symulacji

Próbowałem jeszcze prostszy scenariusz, model regresji przechwycenia tylko dla . Ten sam problem występuje. $Y$

# note: simulation takes 5-10 min on my laptop; can reduce boot.reps
#  and n.sims.run if wanted
# set the number of cores: can change this to match your machine
library(doParallel)
registerDoParallel(cores=8)
boot.reps = 1000
n.sims.run = 1000

for ( j in 1:n.sims.run ) {

  # make initial dataset from which to bootstrap
  # generate under null
  d = data.frame( X1 = rnorm( n = 1000 ), Y1 = rnorm( n = 1000 ) )

  # fit OLS to original data
  mod.orig = lm( Y1 ~ X1, data = d )
  bhat = coef( mod.orig )[["X1"]]
  se = coef(summary(mod.orig))["X1",2]
  rej = coef(summary(mod.orig))["X1",4] < 0.05

  # run all bootstrap iterates
  parallel.time = system.time( {
    r = foreach( icount( boot.reps ), .combine=rbind ) %dopar% {

      # Algorithm 6.2: Resample entire cases - FAILS
      # residuals of this model are repeated, so not normal?
      ids = sample( 1:nrow(d), replace=TRUE )
      b = d[ ids, ]

      # # Method 2: Resample just the residuals themselves - WORKS
      # b = data.frame( X1 = d$X1, Y1 = sample(mod.orig$residuals, replace = TRUE) )

      # # Method 3: Subsampling without replacement - WORKS
      # ids = sample( 1:nrow(d), size = 500, replace=FALSE )
      # b = d[ ids, ]

      # save stats from bootstrap sample
      mod = lm( Y1 ~ X1, data = b ) 
      data.frame( bhat = coef( mod )[["X1"]],
                  se = coef(summary(mod))["X1",2],
                  rej = coef(summary(mod))["X1",4] < 0.05 )

    }
  } )[3]


  ###### Results for This Simulation Rep #####
  r = data.frame(r)
  names(r) = c( "bhat.bt", "se.bt", "rej.bt" )

  # return results of each bootstrap iterate
  new.rows = data.frame( bt.iterate = 1:boot.reps,
                         bhat.bt = r$bhat.bt,
                         se.bt = r$se.bt,
                         rej.bt = r$rej.bt )
  # along with results from original sample
  new.rows$bhat = bhat
  new.rows$se = se
  new.rows$rej = rej

  # add row to output file
  if ( j == 1 ) res = new.rows
  else res = rbind( res, new.rows )
  # res should have boot.reps rows per "j" in the for-loop

  # simulation rep counter
  d$sim.rep = j

}  # end loop over j simulation reps



##### Analyze results #####

# dataset with only one row per simulation
s = res[ res$bt.iterate == 1, ]

# prob of rejecting within each resample
# should be 0.05
mean(res$rej.bt); mean(s$rej)

Aktualizacja nr 2: odpowiedź

W komentarzach i odpowiedziach zaproponowano kilka możliwości i wykonałem więcej symulacji, aby je empirycznie przetestować. Okazuje się, że JWalker ma rację, że problem polega na tym, że musieliśmy wyśrodkować statystyki bootstrapu według szacunków oryginalnych danych, aby uzyskać prawidłowy rozkład próbkowania w . Myślę jednak również, że komentarz Whubera dotyczący naruszenia założeń testu parametrycznego jest również poprawny, chociaż w tym przypadku faktycznie otrzymujemy nominalne fałszywe alarmy, gdy naprawimy problem JWalkera. $H_0$

r bootstrap simulation resampling półprzejście
źródło

1

W standardowym bootstrapie rozważa się jedynie rozkład bootstrapu współczynnika X1, a nie powiązane z nim wartości p. Dlatego nie jest to problem z bootstrapem. Niemniej jednak twoja obserwacja jest interesująca i nieintuicyjna.

Michael M

1

@MichaelM, to prawda. Ale ponieważ wspólny CDF danych w próbkach powinien zbiegać się w n, a liczba bootstrapów iteruje do prawdziwego CDF, który wygenerował oryginalne dane, nie spodziewałbym się, że wartości p będą się również różnić.

półfinał

Dobrze. Jestem całkiem pewien, że efekt pochodzi z obserwacji, które nie są niezależne (jak powiedziałeś), powodując zbyt optymistyczne błędy standardowe. W twojej symulacji wydaje się to jedynym naruszonym założeniem normalnego modelu liniowego. Może uda nam się nawet uzyskać odpowiedni współczynnik deflacji wariancji.

Michael M

2

Jedną z rzeczy, które są jasne w metodzie 1, jest naruszenie założenia błędu iid: podczas ponownego próbkowania z zamianą, resztki dla dowolnego

x

$x$ idsids <- unique(ids)

2

@whuber. Widzę. To by tłumaczyło, dlaczego ponowne próbkowanie reszt za pomocą prac zastępczych pomimo powtarzających się obserwacji: resztki tego modelu są ponownie niezależne od X. Jeśli chcesz zrobić z tego odpowiedź, chętnie to zaakceptuję.

półfinał

5

Po ponownym próbkowaniu wartości null oczekiwana wartość współczynnika regresji wynosi zero. Podczas ponownego próbkowania niektórych obserwowanych zestawów danych oczekiwaną wartością jest obserwowany współczynnik dla tych danych. Nie jest to błąd typu I, jeśli P <= 0,05 podczas ponownego próbkowania obserwowanych danych. W rzeczywistości jest to błąd typu II, jeśli masz P> 0,05.

Można uzyskać intuicję, obliczając korelację między abs (b) a średnią (P). Oto prostszy kod do zreplikowania tego, co zrobiłeś, a także obliczenie korelacji między b i błędem „typu I” w zestawie symulacji

boot.reps = 1000
n.sims.run = 10
n <- 1000
b <- matrix(NA, nrow=boot.reps, ncol=n.sims.run)
p <- matrix(NA, nrow=boot.reps, ncol=n.sims.run)
for(sim_j in 1:n.sims.run){
  x <- rnorm(n)
  y <- rnorm(n)
  inc <- 1:n
  for(boot_i in 1:boot.reps){
    fit <- lm(y[inc] ~ x[inc])
    b[boot_i, sim_j] <- abs(coefficients(summary(fit))['x[inc]', 'Estimate'])
    p[boot_i, sim_j] <- coefficients(summary(fit))['x[inc]', 'Pr(>|t|)']
    inc <- sample(1:n, replace=TRUE)
  }
}
# note this is not really a type I error but whatever
type1 <- apply(p, 2, function(x) sum(x <= 0.05))/boot.reps
# correlation between b and "type I"
cor(b[1, ], type1)

aktualizacja odpowiedzi przez grand_chat nie jest powodem, dla którego częstotliwość P <= 0,05 wynosi> 0,05. Odpowiedź jest bardzo prosta i to, co powiedziałem powyżej - oczekiwana wartość średniej każdej ponownej próby jest oryginalną, zaobserwowaną średnią. Jest to cała podstawa bootstrapu, który został opracowany w celu generowania standardowych błędów / limitów ufności na obserwowanej średniej, a nie jako test hipotez. Ponieważ oczekiwanie nie jest równe zero, oczywiście „błąd typu I” będzie większy niż alfa. I dlatego będzie korelacja między wielkością współczynnika (jak daleko od zera) a wielkością odchylenia „błędu typu I” od alfa.

JWalker
źródło

H_{0} : β = \hat{β}

$H_0: \beta = \widehat{\beta}$

H_{0} : β = 0

$H_0: \beta = 0$

H 0 : β = β ˆ

$H0:β=βˆ$

H 0 : β = 0

$H0:β=0$

H 0 : β = 0

$H0:β=0$ we wczesnych etapach badań w celu odfiltrowania przywry, gdy nie masz wystarczającej wiedzy, aby zdefiniować alternatywną hipotezę, a gdy wiadomo więcej, warto przejść do testowania poprawności swojej wiedzy.

ReneBt

2

Jeśli próbujesz z zamiennikiem oryginalnej oryginalnej próbki, wynikowa próbka ładowania początkowego nie jest normalna . Wspólny rozkład próbki bootstrap przebiega zgodnie z rozkładem mieszanki, który najprawdopodobniej będzie zawierał zduplikowane rekordy, podczas gdy zduplikowane wartości mają prawdopodobieństwo zerowe, gdy pobierzesz próbki z normalnego rozkładu.

Jako prosty przykład, jeśli twoja oryginalna próbka to dwie obserwacje z jednoczynnikowego rozkładu normalnego, wówczas próbka ładowania początkowego z wymianą będzie w połowie składać się z oryginalnej próbki, a w połowie będzie składać się z jednej z oryginalnych wartości, zduplikowanych. Oczywiste jest, że wariancja próbki bootstrap będzie średnio mniejsza niż w oryginale - w rzeczywistości będzie to połowa oryginału.

$p$ $t$

grand_chat
źródło

Ale jeśli tak jest, to czy nie mielibyśmy dokładnie tego samego problemu przy ponownym próbkowaniu reszt z wymianą? Jednak w rzeczywistości takie podejście odrzuca z nominalnym prawdopodobieństwem.

półfinał

Ponadto test t przy n = 1000 nie powinien mieć problemu z danymi niestandardowymi.

półfinał

0

Całkowicie zgadzam się z odpowiedzią @ JWalker.

Jest jeszcze jeden aspekt tego problemu. To jest proces ponownego próbkowania. Można oczekiwać, że współczynnik regresji się wokół zera, ponieważ można zakładać Xi Ysą niezależne. Jednak robisz to podczas ponownego próbkowania

ids = sample( 1:nrow(d), replace=TRUE )
  b = d[ ids, ]

co tworzy korelację, ponieważ próbkujesz Xi Yrazem. Załóżmy na przykład, że pierwszym wierszem zestawu danych djest (x1, y1): W zestawie danych o ponownym próbkowaniu P(Y = y1|X = x1) = 1, a jeśli Xi Ysą niezależne,P(Y|X = x1) ma rozkład normalny.

Więc innym sposobem na naprawę tego jest użycie

b = data.frame( X1 = rnorm( n = 1000 ), Y1 = rnorm( n = 1000 ) )

ten sam kod, którego używasz do generowania d , aby uniezależnić X i Y od siebie.

Ten sam powód wyjaśnia, dlaczego działa z resztkowym próbkowaniem (ponieważ Xjest niezależny od nowego Y).

Tianxia Zhou
źródło

Przez pewien czas myślałem również, że spostrzeżone obserwacje mogą być nie-niezależne, ale po zastanowieniu się nad tym o wiele bardziej tak naprawdę nie sądzę, że tak jest: stats.stackexchange.com/questions/339237/...

połowa -przejść

Problem, który opisuję powyżej, różni się od twojego postu. Mówiłeś o niezależności x's. Odniosłem się do korelacji między Xs i Ys.

Tianxia Zhou,

-1

Największym problemem jest to, że wyniki modelu są fałszywe, a zatem bardzo niestabilne, ponieważ model po prostu dopasowuje hałas. W bardzo dosłownym sensie. Y1 nie jest zmienną zależną ze względu na sposób wygenerowania danych próbki.

Edytuj, w odpowiedzi na komentarze. Spróbujmy jeszcze raz wyjaśnić swoje myślenie.

W przypadku OLS ogólną intencją jest odkrywanie i kwantyfikowanie podstawowych relacji w danych. W przypadku rzeczywistych danych zwykle nie znamy ich dokładnie.

Ale to sztuczna sytuacja testowa. Znamy mechanizm generowania danych EXACT, znajduje się on dokładnie w kodzie opublikowanym przez PO

X1 = rnorm (n = 1000), Y1 = rnorm (n = 1000)

Jeśli wyrażymy to w znanej formie regresji OLS, tj

Y1 = przechwycenie + Beta1 * X1 + Błąd,
który staje się
Y1 = średnia (X1) + 0 (X1) + Błąd

Moim zdaniem jest to model wyrażony w FORMIE liniowej, ale tak naprawdę NIE jest to relacja / model liniowy, ponieważ nie ma nachylenia. Beta1 = 0,000000.

Kiedy wygenerujemy 1000 losowych punktów danych, wykres rozproszenia będzie wyglądał jak klasyczny okrągły spray do strzelby. Może istnieć pewna korelacja między X1 i Y1 w konkretnej próbce 1000 losowych punktów, które zostały wygenerowane, ale jeśli tak, to jest to przypadkowa sytuacja. Jeśli OLS znajdzie korelację, tzn. Odrzuci hipotezę zerową, że nie ma nachylenia, ponieważ ostatecznie wiemy, że tak naprawdę nie ma żadnego związku między tymi dwiema zmiennymi, wówczas OLS dosłownie znalazł wzór w Komponencie błędu. Scharakteryzowałem to jako „dopasowanie hałasu” i „fałszywe”.

Ponadto jednym ze standardowych założeń / wymagań OLS jest to, że (+/-) „model regresji liniowej ma„ parametry liniowe ”. Biorąc pod uwagę dane, uważam, że nie spełniamy tego założenia. Stąd leżące u podstaw statystyki dotyczące istotności są niedokładne. Uważam, że naruszenie założenia liniowości jest bezpośrednią przyczyną nieintuicyjnych wyników bootstrapu.

Kiedy po raz pierwszy przeczytałem ten problem, nie zapadło się w to, że OP zamierzał przetestować pod zerową [hipotezą].

Gdyby jednak zestaw danych został wygenerowany, takie same wyniki nie byłyby intuicyjne

X1 = rnorm (n = 1000), Y1 = X1 * .4 + rnorm (n = 1000)?

Doug Dame
źródło

4

Y_{1}

$Y_1$

Y_{1}

$Y_1$

(-1) z tych samych powodów, które podał @whuber.

półfinał

1

Odpowiedź na ostatnie pytanie w twojej edycji: tak, zdecydowanie. Wypróbuj sam z symulacją. (Ale uważaj na interpretację, ponieważ musisz wziąć pod uwagę, co to jest zero i jaki jest prawdziwy stan rzeczy.)

whuber

Dlaczego testy hipotez na ponownie próbkowanych zestawach danych zbyt często odrzucają wartość zerową?

Aktualizacja nr 1: Więcej symulacji

Aktualizacja nr 2: odpowiedź

Odpowiedzi: