Por que a coleta de dados até obter um resultado significativo aumenta a taxa de erro do tipo I?

60

Eu queria saber exatamente por que a coleta de dados até que um resultado significativo (por exemplo, ) seja obtido (por exemplo, p-hacking) aumenta a taxa de erro do tipo I?p<.05

Eu também apreciaria muito uma Rdemonstração desse fenômeno.

Reza
fonte
6
Você provavelmente quer dizer "p-hacking", porque "harking" refere-se a "Hipotese após resultados conhecidos" e, embora isso possa ser considerado um pecado relacionado, não é sobre o que você está perguntando.
whuber
2
Mais uma vez, o xkcd responde a uma boa pergunta com fotos. xkcd.com/882
Jason
7
@ Jason Eu tenho que discordar do seu link; que não fala sobre coleta cumulativa de dados. O fato de mesmo a coleta cumulativa de dados sobre a mesma coisa e o uso de todos os dados necessários para calcular o valor- estar errado é muito mais trivial do que o caso no xkcd. p
Jik
11
@JiK, chamada justa. Eu estava focado no aspecto "continue tentando até obtermos um resultado que gostamos", mas você está absolutamente correto, há muito mais na questão em questão.
27417 Jason
O @whuber e o user163778 deram respostas muito semelhantes às discutidas no caso praticamente idêntico de "Teste A / B (sequencial)" neste tópico: stats.stackexchange.com/questions/244646/… Ali, discutimos em termos de Family Wise Error taxas e necessidade de ajuste do valor p em testes repetidos. De fato, essa questão pode ser encarada como um problema de teste repetido!
Tomka #

Respostas:

87

O problema é que você está se dando muitas chances de passar no teste. É apenas uma versão chique deste diálogo:

Vou virar você para ver quem paga o jantar.

OK, eu chamo de chefes.

Ratos, você ganhou. Melhores dois dos três?


Para entender isso melhor, considere um modelo simplificado - mas realista - desse procedimento seqüencial . Suponha que você comece com um "teste" de um certo número de observações, mas esteja disposto a continuar experimentando por mais tempo para obter um valor-p menor que . A hipótese nula é que cada observação vem (independentemente) de uma distribuição normal padrão. A alternativa é que o venha independentemente de uma distribuição normal de variação de unidade com uma média diferente de zero. A estatística do teste será a média de todas as observações, , divididas pelo erro padrão, . Para um teste bilateral, os valores críticos são osX i X i n ˉ X 1 / 0.05XiXinX¯ 0,0250,975Zα=±1,961/n0.025 e pontos percentuais da distribuição normal padrão, aproximadamente.0.975Zα=±1.96

Este é um bom teste - para um único experimento com um tamanho fixo de amostra . Ele tem exatamente uma chance de de rejeitar a hipótese nula, não importa qual seja .5 % nn5%n

Vamos converter algebricamente isso para um teste equivalente com base na soma de todos os valores,S N = X 1 + X 2 + + X n = n ˉ X .n

Sn=X1+X2++Xn=nX¯.

Assim, os dados são "significativos" quando

|Zα||X¯1/n|=|Snn/n|=|Sn|/n;

isso é,

(1)|Zα|n|Sn|.

Se formos espertos, reduziremos nossas perdas e desistiremos quando crescer muito e os dados ainda não tiverem entrado na região crítica.n

Isso descreve um passeio aleatório . A fórmula equivale a erguer uma "cerca" ou barreira parabólica curva ao redor do gráfico da caminhada aleatória : o resultado é "significativo" se algum ponto da caminhada aleatória atingir a cerca.Sn(1)(n,Sn)

É uma propriedade de passeios aleatórios que, se esperarmos o suficiente, é muito provável que em algum momento o resultado pareça significativo.

Aqui estão 20 simulações independentes até um limite de amostras. Todos eles começam a testar em amostras, momento em que verificamos se cada ponto está fora das barreiras que foram traçadas de acordo com a fórmula . A partir do momento em que o teste estatístico é o primeiro "significativo", os dados simulados são coloridos em vermelho.n=5000n=30(1)

Figura

Você pode ver o que está acontecendo: a caminhada aleatória sobe e desce cada vez mais à medida que aumenta. As barreiras estão se espalhando na mesma proporção - mas nem sempre são rápidas o suficiente para evitar a caminhada aleatória.n

Em 20% dessas simulações, foi encontrada uma diferença "significativa" - geralmente bem cedo - embora em todas elas a hipótese nula esteja absolutamente correta! A execução de mais simulações desse tipo indica que o tamanho real do teste está próximo de em vez do valor pretendido de : ou seja, sua disposição de continuar procurando "significância" até um tamanho de amostra de oferece chance de rejeitar o nulo, mesmo quando o nulo for verdadeiro.25%α=5%500025%

Observe que nos quatro casos "significativos", à medida que os testes continuavam, os dados pararam de parecer significativos em alguns momentos. Na vida real, um pesquisador que para mais cedo perde a chance de observar essas "reversões". Essa seletividade através da parada opcional influencia os resultados.

Nos testes seqüenciais de honestidade e bondade, as barreiras são linhas. Eles se espalham mais rápido do que as barreiras curvas mostradas aqui.

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)
whuber
fonte
12
+1. Alguma caminhada aleatória acaba atravessando as barreiras com probabilidade 1? Eu sei que a distância esperada após etapas é e procurei agora que a constante de proporcionalidade é , que é menor que 1,96. Mas não sei o que fazer disso. nO(n)2/π
Ameba diz Reinstate Monica
10
@amoeba Essa é uma ótima pergunta, que eu tentei o meu melhor para evitar :-). Se eu pudesse calcular a resposta rapidamente (ou soubesse de imediato), eu a teria postado. Infelizmente, estou muito ocupado para abordar analiticamente agora. A simulação mais longa que fiz foi de 1.000 iterações, observando até com . A proporção de resultados "significativos" parece estabilizar perto de . n=5,000,000α=0.051/4
whuber
4
A questão sobre a probabilidade de atingir o limite é interessante. Imagino que a teoria de Einstein sobre o movimento browniano, relacionando-a a uma equação de difusão, possa ser um ângulo interessante. Temos uma função de distribuição que se espalha com uma taxa e uma "perda de partículas" igual à metade do valor da função de distribuição nesse limite (metade se afasta de zero, além da borda, a outra metade remonta). À medida que essa função de distribuição se espalha e diminui, a "perda" se torna menor. Eu imagino que isso efetivamente criará um limite, isto é, 1/4. α=0.05n
Sextus Empiricus 27/10
6
Razão intuitiva por isso que você vai ter um em algum momento quase certamente: Vamos e . O valor após os primeiros ensaios é praticamente independente do valor após os primeiros ensaios. Então você vai ter um número infinito de "quase" independentes -Valores, assim que um deles é garantido para ser . Obviamente, a convergência real é muito mais rápida do que esse argumento diz. (E se você não gosta de , pode tentar ou ...)n 1 = 10 n k + 1 = 10 n k p n k + 1 p n k p < 0,05 10 n k A ( n k ) B B ( n kp<0.05n1=10nk+1=10nkpnk+1pnkp<0.0510nkA(nk)BB(nk)
JiK
10
@CL. I antecipou sua objeção há vários anos: 17 é a minha semente público. De fato, nos primeiros ensaios (muito mais longos) eu obtive consistentemente maiores taxas de significância substancialmente maiores que 20%. Coloquei a semente em 17 para criar a imagem final e fiquei desapontado por o efeito não ser tão dramático. É a vida. Uma postagem relacionada (ilustrando seu ponto) está em stats.stackexchange.com/a/38067/919 .
whuber
18

As pessoas que são novas no teste de hipóteses tendem a pensar que, uma vez que o valor de p fique abaixo de 0,05, a adição de mais participantes diminuirá ainda mais o valor de p. Mas isso não é verdade. Sob a hipótese nula, o valor de p é distribuído uniformemente entre 0 e 1 e pode saltar bastante nesse intervalo.

Simulei alguns dados em R (minhas habilidades em R são bastante básicas). Nesta simulação, coleciono 5 pontos de dados - cada um com uma associação de grupo selecionada aleatoriamente (0 ou 1) e cada um com uma medida de resultado selecionada aleatoriamente ~ N (0,1). A partir do participante 6, conduzo um teste t a cada iteração.

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
}

Os valores de p estão nesta figura. Observe que encontro resultados significativos quando o tamanho da amostra está entre 70 e 75. Se eu parar por aí, acabarei acreditando que minhas descobertas são significativas, porque perderei o fato de que meus valores de p retornaram com uma amostra maior (isso realmente aconteceu comigo uma vez com dados reais). Como sei que ambas as populações têm uma média de 0, isso deve ser um falso positivo. Esse é o problema da adição de dados até p <0,05. Se você adicionar testes suficientes, p acabará por ultrapassar o limite 0,05 e você poderá encontrar um efeito significativo em qualquer conjunto de dados.

insira a descrição da imagem aqui

TPM
fonte
11
Obrigado, mas seu Rcódigo não é executado.
Reza
3
@Reza você precisaria criar dfprimeiro (de preferência no tamanho final). Como o código começa a ser escrito na linha 6, a implicação (que se encaixa no texto da resposta) é que o df já existe com 5 linhas já preenchidas. Talvez algo parecido com isto tenha sido planejado: n150<-vector("numeric",150); df<-data.frame(gp=n150,val=n150,pval=n150); init<-1:5; df[init,1]<-c(0,1,0,1,0); df[init,2]<-rnorm(5)(execute o código acima) e talvez: plot(df$pv[6:150])
Glen_b 27/10
@ user263778 resposta muito útil e pertinente muito focada. Mas há muita confusão sobre a interpretação do valor-p chamado beleza da dança.
Subhash C. Davar
@ user163778 - você deve incluir o código para inicializar tudo demasiado
Dason
17

Essa resposta diz respeito apenas à probabilidade de obter um resultado "significativo" e à distribuição do tempo para esse evento no modelo do @ whuber.

Como no modelo de @whuber, denota o valor da estatística de teste após a coleta das observações e assume que as observações são seu padrão normal . Então modo que se comporte como um movimento browniano padrão de tempo contínuo, se ignorarmos por enquanto o fato de termos um processo de tempo discreto (gráfico à esquerda abaixo).S(t)=X1+X2++XttX1,X2,

(1)S(t+h)|S(t)=s0N(s0,h),
S(t)

Seja o primeiro tempo de passagem de através das barreiras dependentes do tempo (o número de observações necessárias antes que o teste se torne significativo).TS(t)±zα/2t

Considere o processo transformado obtido escalando por seu desvio padrão no tempo deixando a nova escala de tempo tal que Segue-se de (1) e (2) que é normalmente distribuído com e Y(τ)S(t)tτ=lnt

(2)Y(τ)=S(t(τ))t(τ)=eτ/2S(eτ).
Y(τ+δ)
E(Y(τ+δ)|Y(τ)=y0)=E(e(τ+δ)/2S(eτ+δ)|S(eτ)=y0eτ/2)(3)=y0eδ/2
Var(Y(τ+δ)|Y(τ)=y0)=Var(e(τ+δ)/2S(eτ+δ)|S(eτ)=y0eτ/2)(4)=1eδ,
isto é, é um processo Ornstein-Uhlenbeck (OU) de média zero com uma variação estacionária de 1 e tempo de retorno 2 (gráfico à direita abaixo).Y(τ)

insira a descrição da imagem aqui

Para o modelo transformado, as barreiras tornam-se constantes independentes do tempo iguais a . Sabe-se então ( Nobile et. Al. 1985 ; Ricciardi & Sato, 1988 ) que o primeiro tempo de passagem do processo OU através dessas barreiras é distribuído aproximadamente exponencialmente com algum parâmetro (dependendo das barreiras em ) (estimadas em para abaixo). Há também uma massa de pontos extra em tamanho in . "Rejeição" de±zα/2TY(τ)λ±zα/2λ^=0.125α=0.05ατ=0H0eventualmente acontece com probabilidade 1. Portanto, (o número de observações que precisam ser coletadas antes de obter um resultado "significativo") segue aproximadamente uma distribuição exponencial de log com o valor esperado Assim, tem uma expectativa finita apenas se (por tempo suficiente) grandes níveis de significância ).T=eT

(5)ET1+(1α)0eτλeλτdτ.
Tλ>1α

O exposto acima ignora o fato de que para o modelo real é discreto e que o processo real é discreto - e não em tempo contínuo. Portanto, o modelo acima superestima a probabilidade de a barreira ter sido atravessada (e subestima o ) porque o caminho da amostra em tempo contínuo pode atravessar a barreira apenas temporariamente entre dois pontos de tempo discretos adjacentes e . Mas tais eventos devem ter probabilidade insignificante para grande . E T t t + 1 tTETtt+1t

A figura a seguir mostra uma estimativa de Kaplan-Meier de na escala log-log, juntamente com a curva de sobrevida para a aproximação exponencial em tempo contínuo (linha vermelha).P(T>t)

insira a descrição da imagem aqui

Código 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)
Jarle Tufto
fonte
Obrigado! Você tem alguma referência (padrão) para esses resultados? Por exemplo, por que o processo Y é um Ornstein-Uhlenbeck e onde podemos encontrar o resultado do tempo de passagem indicado?
Grassie 28/10
11
Não vi essa transformação em nenhum outro lugar, mas acredito (3) e (4) que segue facilmente de (1) e (2) e a normalidade caracteriza completamente o processo da OU. O Google scholar retorna muitos resultados sobre a exponencialidade aproximada das distribuições de tempo da primeira passagem para o processo da OU. Mas acredito que neste caso (dentro da aproximação de tempo contínuo) é exatamente exponencialmente distribuído (exceto por uma massa de pontos extra em ) porque vem da distribuição estacionária do processo . τ = 0 Y ( 0 )Tτ=0Y(0)
Jarle Tufto
@Grassie Consulte também math.stackexchange.com/questions/1900304/…
Jarle Tufto
@Grassie Na verdade, meu argumento baseado na falta de memória foi falha. As durações das excursões fora dos limites não são distribuídas exponencialmente. Portanto, com base no mesmo argumento de stats.stackexchange.com/questions/298828/… , mesmo que venha da distribuição estacionária, o tempo da primeira passagem não é exatamente exponencialmente distribuído. Y(0)
Jarle Tufto
5

É preciso dizer que a discussão acima é para uma visão de mundo freqüentista, cuja multiplicidade vem das chances de você dar dados mais extremas, não das chances de dar um efeito à existência. A causa raiz do problema é que os valores de p e os erros do tipo I usam condicionamento de fluxo de informações de retorno ao passado , o que torna importante "como você chegou aqui" e o que poderia ter acontecido. Por outro lado, o paradigma bayesiano codifica ceticismo sobre um efeito no próprio parâmetro, não nos dados. Isso faz com que cada probabilidade posterior seja interpretada da mesma forma, independentemente de você ter calculado outra probabilidade posterior de um efeito há 5 minutos ou não. Mais detalhes e uma simulação simples podem ser encontrados em http://www.fharrell.com/2017/10/continuous-learning-from-data-no.

Frank Harrell
fonte
11
Vamos imaginar um laboratório liderado pelo Dr. B, que é um devoto bayesiano. O laboratório está estudando a preparação social e produziu um fluxo constante de trabalhos mostrando vários efeitos da preparação, sempre apoiados pelo fator Bayes BF> 10. Se eles nunca fazem testes seqüenciais, parece bastante convincente. Mas digamos que eu saiba que eles sempre fazem testes seqüenciais e continuam recebendo novos assuntos até obterem um AM> 10 em favor dos efeitos primários . Então claramente todo esse trabalho é inútil. O fato de fazerem testes sequenciais + seleção faz uma enorme diferença, não importa se é baseado nos valores de p de BF.
Ameba diz Reinstate Monica
11
Não uso os fatores de Bayes. Mas, se tivessem usado probabilidades posteriores e conduzido cada experimento até a probabilidade posterior de um efeito positivo , não haveria absolutamente nada de errado com essas probabilidades. Veja a citação no início do artigo do meu blog - veja o link na minha resposta acima. O grau de crença sobre o efeito de primer vem dos dados e da crença anterior. Se você (como eu) é muito duvidoso desses efeitos primários, é melhor usar um prévio bastante cético ao calcular as probabilidades posteriores. É isso aí. 0.95
Frank Harrell
11
Eu li o seu post no blog, observei a citação e olhei para um artigo semelhante ( Parada opcional: não há problema para os bayesianos ) que alguém ligou nos comentários a outra resposta. Eu ainda não entendi. Se o "nulo" (ausência de efeitos primários) for verdadeiro, se o Dr. B estiver disposto a amostrar por tempo suficiente, ele poderá obter probabilidade posterior> 0,95 toda vez que executar o experimento (exatamente como o Dr. F seria capaz de obter p <0,05 todas as vezes). Se isso é "absolutamente nada de errado", então não sei o que é.
Ameba diz Reinstate Monica
2
Bem, eu discuto esse "ponto maior". Eu não acho que isso seja verdade. Enquanto eu continuo repetindo, sob o efeito nulo de zero e com um dado anterior (digamos, um amplo anterior contínuo centrado em zero), a amostragem repetida sempre mais cedo ou mais tarde rende> 0,98 probabilidade posterior concentrada acima de zero. Uma pessoa que faz amostragem até que isso aconteça (por exemplo, aplicando esta regra de parada), estará errada todas as vezes . Como você pode dizer que essa pessoa estará errada apenas 0,02 das vezes? Eu não entendo Sob essas circunstâncias particulares, não, ele não vai, ele sempre estará errado.
Ameba diz Reinstate Monica
2
Eu não acho que sou. Meu argumento principal é que é injusto e inconsistente culpar simultaneamente os procedimentos freqüentadores por sofrerem testes sequenciais e defender os procedimentos bayesianos como não afetados pelos testes sequenciais. Meu argumento (que é um fato matemático) é que ambos são afetados exatamente da mesma maneira, o que significa que o teste seqüencial pode aumentar o erro tipo I Bayesiano até 100%. Obviamente, se você diz que, por princípio, não está interessado nas taxas de erro do tipo I, isso é irrelevante. Mas os procedimentos freqüentistas também não devem ser responsabilizados por isso.
Ameba diz Reinstate Monica
3

Consideramos um pesquisador coletando uma amostra do tamanho , , para testar alguma hipótese . Ele rejeita se uma estatística de teste adequada exceder seu valor crítico de nível . Caso contrário, ele coletará outra amostra de tamanho , e rejeitará se o teste rejeitar a amostra combinada . Se ele ainda não obtiver rejeição, procederá dessa maneira, até vezes no total.x 1 θ = θ 0 t α c n x 2 ( x 1 , x 2 ) Knx1θ=θ0tαcnx2(x1,x2)K

Este problema parece já ter sido abordado por P. Armitage, CK McPherson e BC Rowe (1969), Jornal da Royal Statistical Society. Série A (132), 2, 235-244: "Testes de significância repetidos em dados acumulativos" .

O ponto de vista bayesiano sobre esta questão, também discutido aqui, é, a propósito, discutido em Berger e Wolpert (1988), "O Princípio da Verossimilhança" , Seção 4.2.

Aqui está uma replicação parcial dos resultados de Armitage et al (código abaixo), que mostra como os níveis de significância aumentam quando , bem como possíveis fatores de correção para restaurar os valores críticos de nível . Observe que a pesquisa na grade demora um pouco para ser executada - a implementação pode ser bastante ineficiente.αK>1α

Tamanho da regra de rejeição padrão em função do número de tentativasK

insira a descrição da imagem aqui

Tamanho em função do aumento de valores críticos para diferentesK

insira a descrição da imagem aqui

Valores críticos ajustados para restaurar 5% dos testes em função deK

insira a descrição da imagem aqui

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")
Christoph Hanck
fonte