Dlaczego gromadzenie danych do momentu uzyskania znaczącego wyniku zwiększa poziom błędu Typu I?

Zastanawiałem się dokładnie, dlaczego gromadzenie danych, dopóki nie zostanie uzyskany znaczący wynik (np. ) (tj. Hakowanie p), zwiększy poziom błędu Typu I?$p \lt .05$

Byłbym również bardzo wdzięczny za Rpokazanie tego zjawiska.

Problem polega na tym, że dajesz sobie zbyt wiele szans na zdanie testu. To tylko fantazyjna wersja tego okna dialogowego:

Odwrócę cię, aby zobaczyć, kto płaci za obiad.

OK, dzwonię po głowy.

Szczury, wygrałeś. Najlepsze dwa z trzech?

---

Aby lepiej to zrozumieć, rozważ uproszczony - ale realistyczny - model tej sekwencyjnej procedury . Załóżmy, że zaczniesz od „serii próbnej” określonej liczby obserwacji, ale jesteś skłonny kontynuować eksperymenty dłużej, aby uzyskać wartość p mniejszą niż . Hipotezą zerową jest to, że każda obserwacja pochodzi (niezależnie) ze standardowego rozkładu normalnego. Alternatywą jest to, że pochodzą niezależnie od rozkładu normalnego wariancji jednostkowej ze średnią niezerową. Statystyka testowa będzie średnią ze wszystkich obserwacji, , podzieloną przez ich błąd standardowy, . W przypadku testu dwustronnego wartościami krytycznymi są:X i X i n ˉ X 1 / $0.05$$X_i$$X_i$$n$$\bar X$ 0,0250,975Zα=±$1/\sqrt{n}$$0.025$ i punktów procentowych standardowego rozkładu normalnego, w przybliżeniu .$0.975$$Z_\alpha=\pm 1.96$

To dobry test - dla pojedynczego eksperymentu ze stałą wielkością próby . Ma dokładnie szans na odrzucenie hipotezy zerowej, bez względu na to, jakie może być .5 % $n$$5\%$$n$

Algebraicznie przekonwertujmy to na równoważny test oparty na sumie wszystkich wartości,S n = X 1 + x 2 + ⋯ + X n = n ˉ X$n$

$$S_n=X_1+X_2+\cdots+X_n = n\bar X.$$

Zatem dane są „znaczące”, kiedy

$$\left| Z_\alpha\right| \le \left| \frac{\bar X}{1/\sqrt{n}} \right| = \left| \frac{S_n}{n/\sqrt{n}} \right| = \left| S_n \right| / \sqrt{n};$$

to jest,

$$\left| Z_\alpha\right| \sqrt{n} \le \left| S_n \right| .\tag{1}$$

Jeśli jesteśmy sprytni, zmniejszymy straty i poddamy się, gdy bardzo duże, a dane nadal nie dotrą do regionu krytycznego.$n$

Opisuje losowy spacer . Wzór sprowadza się do wzniesienia zakrzywionego parabolicznego „ogrodzenia” lub bariery wokół wykresu losowego spaceru : wynik jest „znaczący”, jeśli którykolwiek punkt losowego spaceru uderzy o ogrodzenie.$S_n$$(1)$$(n, S_n)$

Właściwością losowych spacerów jest to, że jeśli będziemy czekać wystarczająco długo, bardzo prawdopodobne jest, że w pewnym momencie wynik będzie wyglądał na znaczący.

Oto 20 niezależnych symulacji do limitu próbek. Wszystkie rozpoczynają testowanie przy próbkach, w tym momencie sprawdzamy, czy każdy punkt znajduje się poza barierami, które zostały narysowane zgodnie ze wzorem . Od momentu, w którym test statystyczny jest po raz pierwszy „znaczący”, symulowane dane są zabarwione na czerwono.$n=5000$$n=30$$(1)$

Możesz zobaczyć, co się dzieje: losowy marsz rośnie i spada w górę i w dół wraz ze wzrostem liczby . Bariery rozprzestrzeniają się w przybliżeniu w tym samym tempie - ale nie wystarczająco szybko, aby uniknąć przypadkowego spaceru.$n$

W 20% tych symulacji wykryto „znaczącą” różnicę - zwykle dość wcześnie - nawet jeśli w każdej z nich hipoteza zerowa jest całkowicie poprawna! Uruchomienie większej liczby symulacji tego typu wskazuje, że rzeczywisty rozmiar testu jest zbliżony do zamiast zamierzonej wartości : oznacza to, że gotowość do szukania „istotności” aż do próby daje szans na odrzucenie wartości null, nawet jeśli wartość null jest prawdziwa.$25\%$$\alpha=5\%$$5000$$25\%$

Zauważ, że we wszystkich czterech „znaczących” przypadkach, gdy testy były kontynuowane, dane przestały wyglądać znaczące w niektórych punktach. W prawdziwym życiu eksperymentator, który zatrzymuje się wcześnie, traci szansę zaobserwowania takich „zmian”. Ta selektywność poprzez opcjonalne zatrzymanie wpływa negatywnie na wyniki.

W sekwencyjnych testach uczciwość do dobroci bariery są liniami. Rozprzestrzeniają się szybciej niż pokazane tutaj zakrzywione bariery.

library(data.table)
library(ggplot2)

alpha <- 0.05   # Test size
n.sim <- 20     # Number of simulated experiments
n.buffer <- 5e3 # Maximum experiment length
i.min <- 30     # Initial number of observations
#
# Generate data.
#
set.seed(17)
X <- data.table(
  n = rep(0:n.buffer, n.sim),
  Iteration = rep(1:n.sim, each=n.buffer+1),
  X = rnorm((1+n.buffer)*n.sim)
)
#
# Perform the testing.
#
Z.alpha <- -qnorm(alpha/2)
X[, Z := Z.alpha * sqrt(n)]
X[, S := c(0, cumsum(X))[-(n.buffer+1)], by=Iteration]
X[, Trigger := abs(S) >= Z & n >= i.min]
X[, Significant := cumsum(Trigger) > 0, by=Iteration]
#
# Plot the results.
#
ggplot(X, aes(n, S, group=Iteration)) +
  geom_path(aes(n,Z)) + geom_path(aes(n,-Z)) +
  geom_point(aes(color=!Significant), size=1/2) +
  facet_wrap(~ Iteration)

Ludzie, którzy są nowi w testowaniu hipotez, mają tendencję do myślenia, że ​​gdy wartość ap spadnie poniżej 0,05, dodanie większej liczby uczestników jeszcze bardziej obniży wartość p. Ale to nie jest prawda. Zgodnie z hipotezą zerową wartość ap jest równomiernie rozłożona między 0 a 1 i może odbijać się dość mocno w tym zakresie.

Symulowałem niektóre dane w R (moje umiejętności R są dość podstawowe). W tej symulacji zbieram 5 punktów danych - każdy z losowo wybranym członkostwem w grupie (0 lub 1) i każdy z losowo wybraną miarą wyniku ~ N (0,1). Począwszy od uczestnika 6, przeprowadzam test t przy każdej iteracji.

for (i in 6:150) {
  df[i,1] = round(runif(1))
  df[i,2] = rnorm(1)
  p = t.test(df[ , 2] ~ df[ , 1], data = df)$p.value
  df[i,3] = p
}

Wartości p są na tym rysunku. Zauważ, że znajduję znaczące wyniki, gdy wielkość próbki wynosi około 70-75. Jeśli się tam zatrzymam, uwierzę, że moje odkrycia są znaczące, ponieważ pominąłem fakt, że moje wartości p podskoczyły z większą próbką (tak naprawdę zdarzyło mi się to raz z prawdziwymi danymi). Ponieważ wiem, że obie populacje mają średnią 0, to musi być fałszywie dodatni. Jest to problem z dodawaniem danych do p <0,05. Jeśli dodasz przeprowadzić wystarczającą liczbę testów, p ostatecznie przekroczy próg 0,05 i możesz znaleźć znaczący efekt dla dowolnego zestawu danych.

Ta odpowiedź dotyczy tylko prawdopodobieństwa uzyskania „istotnego” wyniku i podziału czasu na to zdarzenie w modelu @ Whubera.

Podobnie jak w modelu @whuber, niech oznacza wartość statystyki testowej po zebraniu obserwacji i załóżmy, że obserwacje są w normie normalne . Wtedy tak, że zachowuje się jak standardowy ruch Browna w czasie ciągłym, jeśli na razie zignorujemy fakt, że mamy proces w czasie dyskretnym (lewy wykres poniżej).$S(t)=X_1 + X_2 + \dots + X_t$$t$$X_1,X_2,\dots$

$$S(t+h)|S(t)=s_0 \sim N(s_0, h), \tag{1}$$ $S(t)$

Niech oznacza czas pierwszego przejścia przez bariery zależne od (liczba obserwacji potrzebnych, zanim test stanie się znaczący).$T$$S(t)$$\pm z_{\alpha/2}\sqrt{t}$

Rozważ przetworzony proces uzyskany przez skalowanie przez jego odchylenie standardowe w czasie oraz przez pozostawienie nowej skali czasu takiej, że Z (1) i (2) wynika, że jest zwykle dystrybuowane z i $Y(\tau)$$S(t)$$t$$\tau=\ln t$

$$Y(\tau)=\frac{S(t(\tau))}{\sqrt{t(\tau)}}=e^{-\tau/2}S(e^\tau). \tag{2}$$ $Y(\tau+\delta)$ $$\begin{align} E(Y(\tau+\delta)|Y(\tau)=y_0) &=E(e^{-(\tau+\delta)/2}S(e^{\tau+\delta})|S(e^\tau)=y_0e^{\tau/2}) \\&=y_0e^{-\delta/2} \tag{3} \end{align}$$ $$\begin{align} \operatorname{Var}(Y(\tau+\delta)|Y(\tau)=y_0) &=\operatorname{Var}(e^{(\tau+\delta)/2}S(e^{\tau+\delta})|S(e^\tau)=y_0e^{\tau/2}) \\&=1-e^{-\delta}, \tag{4} \end{align}$$ to znaczy, jest zerowym średnim procesem Ornsteina-Uhlenbecka (OU) ze stacjonarną wariancją 1 i czasem powrotu 2 (prawy wykres poniżej).$Y(\tau)$

W przypadku modelu przekształconego bariery stają się stałymi niezależnymi od czasu równymi . Wiadomo następnie ( Nobile i in. 1985 ; Ricciardi i Sato, 1988 ), że pierwszy czas przejścia procesu OU przez te bariery jest w przybliżeniu wykładniczo rozkładany z pewnym parametrem (w zależności od barier w ) (szacowany na dla poniżej). Istnieje również dodatkowa masa punktowa o rozmiarze w . „Odrzucenie”$\pm z_{\alpha/2}$$\mathcal{T}$$Y(\tau)$$\lambda$$\pm z_{\alpha/2}$$\hat\lambda=0.125$$\alpha=0.05$$\alpha$$\tau=0$$H_0$ostatecznie dzieje się z prawdopodobieństwem 1. Stąd (liczba obserwacji, które należy zebrać przed uzyskaniem „znaczącego” wyniku) w przybliżeniu odpowiada logarytmicznemu rozkładowi wykładniczemu o oczekiwanej wartości Zatem ma skończone oczekiwanie tylko wtedy, gdy (wystarczająco duże poziomy znaczenia ).$T=e^\mathcal{T}$

$$ET\approx 1+(1-\alpha)\int_0^\infty e^\tau \lambda e^{-\lambda \tau}d\tau.\tag{5}$$ $T$$\lambda>1$$\alpha$

Powyższe ignoruje fakt, że dla modelu rzeczywistego jest dyskretny, a prawdziwy proces jest dyskretny, a nie ciągły. Dlatego powyższy model przecenia prawdopodobieństwo przekroczenia bariery (i nie docenia ), ponieważ ścieżka próbki w czasie ciągłym może przekroczyć barierę tylko tymczasowo pomiędzy dwoma sąsiadującymi dyskretnymi punktami czasowymi i . Ale takie zdarzenia powinny mieć znikome prawdopodobieństwo dla dużego . E T t t + 1 $T$$ET$$t$$t+1$$t$

Poniższy rysunek pokazuje oszacowanie Kaplana-Meiera na skali log-log wraz z krzywą przeżycia dla wykładniczego przybliżenia ciągłego czasu (czerwona linia).$P(T>t)$

Kod R:

# Fig 1
par(mfrow=c(1,2),mar=c(4,4,.5,.5))
set.seed(16)
n <- 20
npoints <- n*100 + 1
t <- seq(1,n,len=npoints)
subset <- 1:n*100-99
deltat <- c(1,diff(t))
z <- qnorm(.975)
s <- cumsum(rnorm(npoints,sd=sqrt(deltat)))
plot(t,s,type="l",ylim=c(-1,1)*z*sqrt(n),ylab="S(t)",col="grey")
points(t[subset],s[subset],pch="+")
curve(sqrt(t)*z,xname="t",add=TRUE)
curve(-sqrt(t)*z,xname="t",add=TRUE)
tau <- log(t)
y <- s/sqrt(t)
plot(tau,y,type="l",ylim=c(-2.5,2.5),col="grey",xlab=expression(tau),ylab=expression(Y(tau)))
points(tau[subset],y[subset],pch="+")
abline(h=c(-z,z))

# Fig 2
nmax <- 1e+3
nsim <- 1e+5
alpha <- .05
t <- numeric(nsim)
n <- 1:nmax
for (i in 1:nsim) {
  s <- cumsum(rnorm(nmax))
  t[i] <- which(abs(s) > qnorm(1-alpha/2)*sqrt(n))[1]
}
delta <- ifelse(is.na(t),0,1)
t[delta==0] <- nmax + 1
library(survival)
par(mfrow=c(1,1),mar=c(4,4,.5,.5))
plot(survfit(Surv(t,delta)~1),log="xy",xlab="t",ylab="P(T>t)",conf.int=FALSE)
curve((1-alpha)*exp(-.125*(log(x))),add=TRUE,col="red",from=1,to=nmax)

Trzeba powiedzieć, że powyższa dyskusja jest dla częstokroć postrzegającego świat, dla którego wielość pochodzi z szans, że dane są bardziej ekstremalne, a nie z szans, że dasz efekt, by istnieć. Główną przyczyną problemu jest to, że wartości p i błędy typu I wykorzystują warunkowanie przepływu informacji wstecz w czasie wstecznym, co sprawia, że ​​ważne jest „jak się tu dostałeś” i co mogło się stać zamiast tego. Z drugiej strony paradygmat bayesowski koduje sceptycyzm co do wpływu na sam parametr, a nie na dane. To sprawia, że ​​każde tylne prawdopodobieństwo należy interpretować tak samo, niezależnie od tego, czy obliczono kolejne tylne prawdopodobieństwo efektu 5 minut temu, czy nie. Więcej informacji i prosta symulacja można znaleźć na stronie http://www.fharrell.com/2017/10/continuous-learning-from-data-no.

Rozważamy badacza pobierającego próbkę o rozmiarze , , aby przetestować hipotezę . Odrzuca, jeśli odpowiednia statystyka testowa przekracza jego wartość krytyczną na poziomie . Jeśli nie, pobiera kolejną próbkę o rozmiarze , i odrzuca, jeśli test odrzuca próbkę łączoną . Jeśli nadal nie otrzyma odrzucenia, postępuje w ten sposób, łącznie do razy.x 1 θ = θ 0 t α c n x 2 ( x 1 , x 2 ) $n$$x_1$$\theta=\theta_0$$t$$\alpha$$c$$n$$x_2$$(x_1,x_2)$$K$

Wydaje się, że ten problem został już rozwiązany przez P. Armitage, CK McPherson i BC Rowe (1969), Journal of the Royal Statistics Society. Seria A (132), 2, 235–244: „Powtarzane testy istotności gromadzonych danych” .

Nawiasem mówiąc omawiany tutaj bayesowski punkt widzenia na ten temat omówiono w Berger i Wolpert (1988), „Zasada wiarygodności” , sekcja 4.2.

Oto częściowa replikacja wyników Armitage'a i in. (Kod poniżej), który pokazuje, jak poziomy istotności zwiększają się, gdy , a także możliwe współczynniki korekcyjne w celu przywrócenia wartości krytycznych poziom . Pamiętaj, że wyszukiwanie siatki zajmuje trochę czasu - implementacja może być raczej nieefektywna.$K>1$$\alpha$

Rozmiar standardowej reguły odrzucania jako funkcja liczby prób$K$

Rozmiar jako funkcja rosnących wartości krytycznych dla różnych$K$

Skorygowano wartości krytyczne, aby przywrócić 5% testów w funkcji$K$

reps <- 50000

K <- c(1:5, seq(10,50,5), seq(60,100,10)) # the number of attempts a researcher gives herself
alpha <- 0.05
cv <- qnorm(1-alpha/2)

grid.scale.cv <- cv*seq(1,1.5,by=.01) # scaled critical values over which we check rejection rates
max.g <- length(grid.scale.cv)
results <- matrix(NA, nrow = length(K), ncol=max.g)

for (kk in 1:length(K)){
  g <- 1
  dev <- 0
  K.act <- K[kk]
  while (dev > -0.01 & g <= max.g){
    rej <- rep(NA,reps)
    for (i in 1:reps){
      k <- 1
      accept <- 1
      x <- rnorm(K.act)
      while(k <= K.act & accept==1){
        # each of our test statistics for "samples" of size n are N(0,1) under H0, so just scaling their sum by sqrt(k) gives another N(0,1) test statistic
        rej[i] <- abs(1/sqrt(k)*sum(x[1:k])) > grid.scale.cv[g] 
        accept <- accept - rej[i]
        k <- k+1
      }
    }
    rej.rate <- mean(rej)
    dev <- rej.rate-alpha
    results[kk,g] <- rej.rate
    g <- g+1
  }
}
plot(K,results[,1], type="l")
matplot(grid.scale.cv,t(results), type="l")
abline(h=0.05)

cv.a <- data.frame(K,adjusted.cv=grid.scale.cv[apply(abs(results-alpha),1,which.min)])
plot(K,cv.a$adjusted.cv, type="l")