Pourquoi la collecte de données jusqu'à l'obtention d'un résultat significatif augmente-t-elle le taux d'erreur de type I?

60

Je me demandais exactement pourquoi la collecte de données jusqu'à l' obtention d' un résultat significatif (par exemple, ) (par exemple, p-piratage) augmente le taux d'erreur de type I?p<.05

J'apprécierais aussi beaucoup une Rdémonstration de ce phénomène.

Reza
la source
6
Vous voulez probablement dire «p-hacking», parce que «harking» fait référence à «hypothéquer après que les résultats sont connus» et, bien que cela puisse être considéré comme un péché connexe, ce n'est pas ce que vous semblez demander.
whuber
2
Encore une fois, xkcd répond à une bonne question par des images. xkcd.com/882
Jason
7
@ Jason, je suis en désaccord avec votre lien; cela ne parle pas de la collecte cumulative de données. Le fait que même la collecte cumulative de données sur la même chose et l' utilisation de toutes les données pour calculer la valeur soit erronée est beaucoup plus anodin que dans le cas de ce xkcd. p
JiK
1
@JiK, appel juste. J'étais concentré sur l'aspect "continue d'essayer jusqu'à ce que nous obtenions un résultat qui nous plait", mais tu as tout à fait raison, il y a beaucoup plus que cela dans la question à l'examen.
Jason
@whuber et user163778 a donné des réponses très similaire à celle décrite pour le cas pratiquement identique de « A / B (séquentielle) test » dans ce fil: stats.stackexchange.com/questions/244646/... Là, nous avons discuté en termes de Wise Family Error les taux et la nécessité d'ajuster la valeur p lors d'essais répétés. En fait, cette question peut être considérée comme un problème de test répété!
tomka

Réponses:

87

Le problème est que vous vous donnez trop de chances de réussir le test. C'est juste une version sophistiquée de cette boîte de dialogue:

Je vais vous retourner pour voir qui paie pour le dîner.

OK, j'appelle des têtes.

Les rats, vous avez gagné. Deux meilleurs sur trois?


Pour mieux comprendre cela, considérons un modèle simplifié - mais réaliste - de cette procédure séquentielle . Supposons que vous commenciez par un "essai" d'un certain nombre d'observations, mais que vous souhaitiez continuer à expérimenter plus longtemps afin d'obtenir une valeur p inférieure à . L'hypothèse nulle est que chaque observation provient (indépendamment) d'une distribution normale standard. L’alternative est que les proviennent indépendamment d’une distribution normale à variance unitaire avec une moyenne non nulle. La statistique de test sera la moyenne de toutes les observations, , divisée par leur erreur standard, . Pour un test bilatéral, les valeurs critiques sont les0.05XiXinX¯1/n0.025 et point de pourcentage de la distribution normale standard, environ.0.975Zα=±1.96

C'est un bon test - pour une seule expérience avec une taille d'échantillon fixe . Il a exactement chances de rejeter l'hypothèse nulle, quel que soit .n5%n

Convertissons-le algébriquement en un test équivalent basé sur la somme de toutes les valeurs,n

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

Ainsi, les données sont "significatives" lorsque

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

C'est,

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

Si nous sommes intelligents, nous réduirons nos pertes et abandonnerons dès que grandira très fort et les données ne sont toujours pas entrées dans la région critique.n

Ceci décrit une marche aléatoire . La formule revient à ériger une "barrière" parabolique incurvée autour du tracé de la marche aléatoire : le résultat est "significatif" si un point quelconque de la marche aléatoire heurte la clôture.Sn(1)(n,Sn)

C’est une propriété des promenades aléatoires que si nous attendons assez longtemps, il est très probable qu’à un moment donné, le résultat paraisse important.

Voici 20 simulations indépendantes jusqu’à une limite de échantillons. Ils commencent tous à tester échantillons. Nous vérifions ensuite si chaque point se situe en dehors des barrières tracées selon la formule . À partir du moment où le test statistique est «significatif», les données simulées sont colorées en rouge.n=5000n=30(1)

Figure

Vous pouvez voir ce qui se passe: la marche aléatoire monte et descend de plus en plus à mesure que augmente. Les barrières se dispersent à peu près au même rythme - mais pas assez vite pour éviter toujours la marche au hasard.n

Dans 20% de ces simulations, une différence «significative» a été observée - généralement assez tôt - même si l'hypothèse nulle est absolument correcte dans chacune d'entre elles! Effectuer plus de simulations de ce type indique que la taille réelle du test est proche de plutôt que la valeur souhaitée de : c'est-à-dire votre volonté de continuer à rechercher une "importance" jusqu'à une taille de vous donne chances de rejeter la valeur null, même si celle-ci est vraie.25%α=5%500025%

Notez que dans les quatre cas « importants », comme poursuite des essais, les données ont arrêté la recherche importante sur certains points. Dans la vie réelle, un expérimentateur qui s’arrête de bonne heure perd la chance d’observer de tels «renversements». Cette sélectivité par arrêt optionnel biaise les résultats.

Dans les tests séquentiels honnêtes, les barrières sont des lignes. Ils se propagent plus rapidement que les barrières courbes montrées ici.

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
la source
12
+1 Est-ce qu'une marche aléatoire donnée traverse éventuellement les barrières avec une probabilité de 1? Je sais que la distance attendue après étapes est et je regarde maintenant que la constante de proportionnalité est , inférieure à 1,96. Mais je ne sais pas quoi en faire. nO(n)2/π
amibe dit de réintégrer Monica
10
@ amoeba C'est une excellente question, que j'ai essayé de mon mieux d'esquiver :-). Si j’avais pu calculer rapidement la réponse (ou si je l’avais su à l’avance), je l’aurais postée. Malheureusement, je suis trop occupé pour y répondre de manière analytique pour le moment. La plus longue simulation que j'ai faite était 1 000 itérations allant jusqu'à avec . La proportion de résultats "significatifs" semble se stabiliser vers . n=5,000,000α=0.051/4
whuber
4
La question sur la probabilité d'atteindre la limite est intéressante. J'imagine que la théorie des Einstein sur le mouvement brownien, la reliant à une équation de diffusion, pourrait constituer un angle intéressant. Nous avons une fonction de distribution s'étalant avec un taux et une "perte de particules" égale à la moitié de la valeur de la fonction de distribution à cette limite (la moitié s'éloigne de zéro, au-delà de la frontière, l'autre moitié il revient). Au fur et à mesure que cette fonction de distribution se répand et s'amincit, la "perte" diminue. J'imagine que cela va effectivement créer une limite, c'est-à-dire ce 1/4. α=0.05n
Sextus Empiricus
6
La raison intuitive pour laquelle vous obtiendrez sûrement un à un moment donné est: Soient et . La valeur après les premiers essais est à peu près indépendante de la valeur après les premiers essais . Donc, vous aurez un nombre infini de valeurs "presque" indépendantes , donc il est garanti que l’une d’elles est . Bien entendu, la convergence réelle est beaucoup plus rapide que ne le dit cet argument. (Et si vous n'aimez pas , vous pouvez essayer ou ...)p<0.05n1=10nk+1=10nkpnk+1pnkp<0.0510nkA(nk)BB(nk)
JiK
10
@CL. J'ai anticipé votre objection il y a plusieurs années: 17 est ma semence publique. En fait, au début des essais (beaucoup plus longs), mon taux de signification était constamment supérieur à 20%. J'ai mis la graine à 17 pour créer l'image finale et j'ai été déçu que l'effet n'ait pas été aussi spectaculaire. C'est la vie. Un article connexe (illustrant votre propos) est à l' adresse stats.stackexchange.com/a/38067/919 .
whuber
18

Les personnes qui débutent dans le test d'hypothèses ont tendance à penser qu'une fois qu'une valeur ap descend en dessous de 0,05, l'ajout de participants ne fera que réduire davantage la valeur p. Mais ce n'est pas vrai. Sous l'hypothèse nulle, la valeur ap est uniformément répartie entre 0 et 1 et peut rebondir assez fort dans cette plage.

J'ai simulé des données dans R (mes compétences en R sont assez basiques). Dans cette simulation, je collecte 5 points de données - chacun avec une appartenance à un groupe sélectionné au hasard (0 ou 1) et chacun avec une mesure de résultat sélectionnée au hasard ~ N (0,1). À partir du participant 6, je réalise un test t à chaque itération.

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
}

Les valeurs p sont dans cette figure. Notez que je trouve des résultats significatifs lorsque la taille de l'échantillon se situe autour de 70-75. Si je m'arrête là, je finirai par croire que mes découvertes sont significatives parce que j'aurai oublié le fait que mes valeurs de p ont remonté avec un échantillon plus grand (cela m'est arrivé en réalité une fois avec des données réelles). Comme je sais que les deux populations ont une moyenne de 0, cela doit être un faux positif. C'est le problème avec l'ajout de données jusqu'à ce que p <0,05. Si vous ajoutez suffisamment de tests, p finira par dépasser le seuil de 0,05 et vous pourrez trouver un effet significatif sur tout jeu de données.

entrez la description de l'image ici

TPM
la source
1
Merci, mais votre Rcode ne fonctionne pas du tout.
Reza
3
@Reza vous devez d' dfabord créer (de préférence à sa taille finale). Depuis le début de l' écriture de code à la ligne 6 l'implication ( ce qui correspond avec le texte de la réponse) est que df existe déjà avec 5 lignes déjà remplies Peut - être quelque chose comme cela était prévu. 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)(Puis exécutez le code ci - dessus) alors peut - être: plot(df$pv[6:150])
Glen_b
@ user263778 très concentré réponse utile et pertinente. Mais il y a trop de confusion sur l'interprétation de la p-valeur ainsi appelée - beauté dansante.
Subhash C. Davar
@ user163778 - vous devriez inclure le code pour tout initialiser aussi
Dason
17

Cette réponse ne concerne que la probabilité d'obtenir finalement un résultat "significatif" et la distribution du temps écoulé avant cet événement selon le modèle de @ whuber.

Comme dans le modèle de @whuber, supposons que désigne la valeur de la statistique de test après que observations aient été collectées et supposons que les observations sont iid standard normal . Alors tels que se comporte comme un mouvement brownien standard à temps continu, si on ignore pour le moment le fait que nous ayons un processus à temps discret (graphique de gauche ci-dessous).S(t)=X1+X2++XttX1,X2,

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

Soit le premier temps de passage de entre les barrières dépendantes du temps (nombre d'observations nécessaires avant que le test ne devienne significatif).TS(t)±zα/2t

Considérons le processus transformé obtenu en mettant à l'échelle par son écart type au temps et en laissant la nouvelle échelle de temps telle que Il résulte de (1) et (2) que est normalement distribué avec et 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δ,
c'est-à-dire que est un processus Ornstein-Uhlenbeck (OU) à moyenne nulle avec une variance stationnaire de 1 et un temps de retour 2 (graphique de droite ci-dessous).Y(τ)

entrez la description de l'image ici

Pour le modèle transformé, les barrières deviennent des constantes indépendantes du temps égales à . On sait alors ( Nobile et al., 1985 ; Ricciardi et Sato, 1988 ) que le premier temps de passage du processus OU OU travers ces barrières est distribué de façon approximativement exponentielle avec un paramètre (en fonction des barrières à ) (estimé à pour ci-dessous). Il existe également une masse de points supplémentaire de taille dans . "Rejet" de±zα/2TY(τ)λ±zα/2λ^=0.125α=0.05ατ=0H0se produit finalement avec la probabilité 1. Par conséquent, (nombre d'observations à collecter avant d'obtenir un résultat "significatif") suit approximativement une distribution log exponentielle avec la valeur attendue Ainsi, n’a une espérance finie que si (suffisamment niveaux importants de signification ).T=eT

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

Ce qui précède ignore le fait que pour le modèle réel est discret et que le processus réel est discret plutôt que continu. Par conséquent, le modèle ci-dessus surestime la probabilité que la barrière soit franchie (et sous-estime ) car le trajet d'échantillonnage en temps continu ne peut franchir la barrière que temporairement entre deux points temporels discrets adjacents et . Mais de tels événements devraient avoir une probabilité négligeable pour les grands . TETtt+1t

La figure suivante montre une estimation de Kaplan-Meier de à l'échelle log-log ainsi que la courbe de survie pour l'approximation à temps continu exponentielle (ligne rouge).P(T>t)

entrez la description de l'image ici

Code 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
la source
Merci! Avez-vous des références (standard) pour ces résultats? Par exemple, pourquoi le processus Y est-il un processus Ornstein-Uhlenbeck et où pouvons-nous trouver le résultat du temps de passage indiqué?
Grassie
1
Je n'ai vu cette transformation nulle part ailleurs, mais je crois (3) et (4) que cela découle facilement de (1) et (2) et que la normalité caractérise complètement le processus de l'OU. Google scholar renvoie beaucoup de résultats sur l'exponentialité approximative des distributions temporelles de premier passage pour le processus d'UO. Mais je crois que dans ce cas (dans l'approximation en temps continu) est exactement distribuée de façon exponentielle (sauf pour une masse de points supplémentaire dans ) car provient de la distribution stationnaire du processus . Tτ=0Y(0)
Jarle Tufto
@Grassie Voir aussi math.stackexchange.com/questions/1900304/...
Jarle Tufto
@ Grassie En fait, mon argument fondé sur l'absence de mémoire était imparfait. Les durées des excursions hors des limites ne sont pas distribuées de manière exponentielle. Par conséquent, basé sur le même argument que dans stats.stackexchange.com/questions/298828/… , même si provient de la distribution stationnaire, le premier temps de passage n’est pas exactement distribué de manière exponentielle. Y(0)
Jarle Tufto
5

Il faut dire que la discussion ci-dessus est pour une vision du monde fréquentiste pour laquelle la multiplicité provient des chances que vous donnez des données d'être plus extrêmes, pas des chances que vous donnez un effet à exister. La cause fondamentale du problème est que les valeurs p et les erreurs de type I utilisent le conditionnement du flux d'informations en amont et en arrière, ce qui rend important "comment vous êtes arrivé ici" et ce qui aurait pu se produire à la place. D'autre part, le paradigme bayésien code le scepticisme quant à un effet sur le paramètre lui-même, pas sur les données. Cela fait que chaque probabilité postérieure soit interprétée de la même manière, que vous calculiez une autre probabilité postérieure d'effet il y a 5 minutes ou non. Plus de détails et une simulation simple peuvent être trouvés à http://www.fharrell.com/2017/10/continuous-learning-from-data-no.

Frank Harrell
la source
1
Imaginons un laboratoire dirigé par le Dr B, un Bayésien dévot. Le laboratoire étudie l’amorçage social et a produit un flux constant d’articles montrant divers effets de l’amorçage, soutenus à chaque fois par le facteur de Bayes BF> 10. S'ils ne font jamais de tests séquentiels, cela paraît assez convaincant. Mais disons que j’apprends qu’ils font toujours des tests séquentiels et continuent à avoir de nouveaux sujets jusqu’à ce qu’ils obtiennent un BF> 10 en faveur des effets d’amorçage . Alors clairement, tout ce travail ne vaut rien. Le fait qu'ils effectuent des tests et une sélection séquentiels fait une différence énorme, peu importe que ce soit basé sur les valeurs p de BF.
amibe dit de réintégrer Monica
1
Je n'utilise pas les facteurs de Bayes. Mais s’ils avaient utilisé des probabilités postérieures et mené chaque expérience jusqu’à la probabilité postérieure d’effet positif , il n’y aurait absolument aucun problème avec ces probabilités. Regardez la citation au début de mon article de blog - voir le lien dans ma réponse ci-dessus. Le degré de conviction concernant l’effet d’amorçage provient des données et de la croyance antérieure. Si vous (comme moi) êtes très sceptique quant à de tels effets d’amorçage, vous feriez mieux d’utiliser un préalable assez sceptique lors du calcul des probabilités postérieures. C'est ça. 0.95
Frank Harrell
1
J'ai lu votre billet de blog, remarqué la citation et consulté un article similaire ( Interruption facultative: aucun problème pour les Bayésiens ) auquel quelqu'un d'autre a associé les commentaires dans une autre réponse. Je ne comprends toujours pas. Si "nul" (effets d'amorçage absents) est vrai, si le Dr B souhaite échantillonner suffisamment longtemps, il sera en mesure d'obtenir une probabilité postérieure> 0,95 à chaque fois qu'il effectue l'expérience (exactement comme le Dr F pourrait le faire). obtenir p <0,05 à chaque fois). Si c'est "absolument rien à redire" alors je ne sais pas ce que c'est.
amibe dit de réintégrer Monica
2
Eh bien, je conteste ce "plus gros point". Je ne pense pas que c'est vrai. Comme je le répète sans cesse, sous l'effet nul de zéro et avec un préalable donné (disons un large préalable continu centré à zéro), un échantillonnage répété donnera toujours tôt ou tard une probabilité> 0.98 de probabilité postérieure concentrée au-dessus de zéro. Une personne qui prélève des échantillons jusqu'à ce que cela se produise (c'est-à-dire en appliquant cette règle d'arrêt) se trompe à chaque fois . Comment pouvez-vous dire que cette personne aura tort que 0,02 du temps? Je ne comprends pas. Dans ces circonstances particulières, non il ne le fera pas, il aura toujours tort.
amibe dit de réintégrer Monica
2
Je ne pense pas que je suis. Mon point le plus important est qu’il est injuste et incohérent d’accuser simultanément les procédures fréquentistes de subir des tests séquentiels et de défendre les procédures bayésiennes comme étant non affectées par les tests séquentiels. Mon point de vue (qui est un fait mathématique) est qu’ils sont tous deux affectés exactement de la même manière, ce qui signifie que des tests séquentiels peuvent augmenter l’erreur bayésienne de type I jusqu’à 100%. Bien sûr, si vous dites que, par principe, les taux d'erreur de type I ne vous intéressent pas, alors c'est sans importance. Mais alors les procédures fréquentistes ne doivent pas être blâmées pour cela non plus.
amibe dit de réintégrer Monica
3

Nous considérons un chercheur rassemblant un échantillon de taille , , pour tester une hypothèse . Il rejette si une statistique de test appropriée dépasse son niveau- valeur critique . Si ce n'est pas le cas, il collecte un autre échantillon de taille , , et le rejette si le test échoue pour l'échantillon combiné . S'il n'obtient toujours aucun rejet, il procède de la sorte, jusqu'à fois au total.nx1θ=θ0tαcnx2(x1,x2)K

P. Armitage, CK McPherson et BC Rowe (1969), Journal de la Royal Statistical Society, ont déjà abordé ce problème . Série A (132), 2, 235-244: "Essais de répétition sur des données accumulées" .

Le point de vue bayésien sur cette question, également discuté ici, est d'ailleurs discuté dans Berger et Wolpert (1988), "Le principe de vraisemblance" , section 4.2.

Voici une réplication partielle des résultats d'Armitage et al. (Code ci-dessous), qui montre comment les niveaux d'importance grossissent lorsque , ainsi que les facteurs de correction possibles pour restaurer les valeurs critiques de niveau . Notez que la recherche dans la grille prend un certain temps - la mise en œuvre peut être plutôt inefficace.K>1α

Taille de la règle de rejet standard en fonction du nombre de tentativesK

entrez la description de l'image ici

Taille en fonction de l'augmentation des valeurs critiques pour différentsK

entrez la description de l'image ici

Valeurs critiques ajustées pour restaurer 5% des tests en fonction deK

entrez la description de l'image ici

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
la source