Motivația
Cu câteva luni în urmă, într-o încăierare secundară în timpul marii controverse privind curba p, Richard McElreath a menționat că valorile p din ipoteza nulă nu sunt întotdeauna distribuite uniform, așa cum se pretinde uneori. Acest lucru m-a determinat să verific fenomenul. Recunosc că am avut în cap ideea de bază că valorile p sunt într-adevăr distribuite uniform dacă ipoteza nulă este adevărată. Se pare că acesta este doar „deseori”, nu întotdeauna.
Așa cum se întâmplă în mod obișnuit pentru acest blog, principala motivație a fost să mă asigur că am înțeles ceva eu însumi, așa că nu există nimic deosebit de nou pentru lume în acest blog. Dar ar putea fi interesant pentru unii. Au fost câteva cârlige ciudate.
Am vrut să verific cazul în care comparăm numărul de „reușite” dintr-un eșantion relativ mic de rezultate binare succes/eșec, împărțit în două grupuri. Ipoteza nulă este că probabilitatea de succes este aceeași în fiecare grup. Am vrut să văd distribuția valorii p pentru un test al acestei ipoteze nule cu dimensiuni ale eșantionului de 10, 100 sau 1000 de observații în fiecare grup; pentru diferite valori ale probabilității de succes subiacente și pentru când dimensiunea eșantionului grupurilor este aleatorie, dar are o medie de 100.
Calcularea valorii p
Cârligele au fost în modul de calcul al valorii p. De fapt, am mers în trei moduri diferite:
- Modul meu de aproximare leneș este să estimez varianța diferenței dintre cele două proporții ale eșantionului în cadrul ipotezei nule și să mă bazez pe normalitatea asimptotică pentru a obține o probabilitate ca valoarea să fie la fel de extremă cum este de fapt. Acest lucru are dezavantajul că se știe că nu este deosebit de corect, mai ales atunci când eșantionul este mic sau probabilitatea de succes este aproape de 1 sau 0. Avantajul este că nu a trebuit să caut nimic și a fost ușor de vectorizat. Această metodă se numește
pval_handîn complotul de mai jos. - O metodă mai bună este să folosești testul exact Fisher, conform celebrei analize lady-tasting-ceai, dar metoda out of the box pe care o foloseam (de la
corporaPachetul R de Stephanie Evert) nu funcționează atunci când succesele observate în ambele eșantioane sunt zero. Această metodă se numeștepval_fisherîn complotul de mai jos. - Cea mai ieșită din cutie dintre toate este utilizarea
prop.test()funcția de lastatspachet integrat în R. Dezavantajul aici a fost deranjarea vectorizării acestei funcții pentru a rula eficient cu un număr mare de simulări. Această metodă se numeștepval_proptestîn complotul de mai jos.
Deoarece eram nervos dacă aceste metode ar putea da rezultate material diferite, am început prin a compara rezultatele pe care le-au dat, cu dimensiuni diferite ale eșantionului și parametri de bază, cu doar 1.000 de repetări pentru fiecare combinație de dimensiune a eșantionului și parametru. Aceasta oferă această comparație:
Deci putem vedea că pval_fisher şi pval_proptest metodele dau efectiv aceleași rezultate, în timp ce metoda mea făcută manual are un număr semnificativ de discrepanțe. Din această cauză am decis să rămân cu Evert’s corpora::fisher.pval. Tocmai l-am întărit cu o funcție wrapper pentru a defini valoarea p (șansa de a vedea date la fel de extreme ca aceasta, dacă ipoteza nulă este adevărată) să fie 1 dacă succesele observate sunt 0 în ambele eșantioane:
library(tidyverse)
library(corpora)
library(GGally)
library(glue)
library(scales)
library(frs) # for svg_png()
#' Version of fisher.pval that won't break if k1 and k2 both 0
tough_fisher <- function(k1, n1, k2, n2, set_both_zero = 1, ...){
problems <- k1 + k2 == 0
k1(problems) <- 1
temp <- corpora::fisher.pval(k1, n1, k2, n2)
temp(problems) <- set_both_zero
return(temp)
}
# note needed a decision on what to do when k1 and k2 are both zero. I think
# p value is 1 here, as p is the "probability of seeing data at least as extreme
# as this if the null hypothesis of no difference is true". Not sure why Fisher's
# exact test returns an error, worth looking into that.
Rezultate
Așa că, cu această problemă din drum, mi-am propus să calculez o mulțime de valori p. Am făcut asta în situații
- unde cele două grupuri au fost de dimensiuni egale cu 10, 100 sau 1.000 de observații fiecare; sau unde sunt un număr aleatoriu de observații, medie 100, tot aceleași pentru fiecare dintre cele două grupuri; şi
- cu o probabilitate de succes de 0,2, 0,5 sau 0,8.
Ideea pentru unele dintre motivele pentru care valorile p nu sunt distribuite uniform este că, cu un eșantion finit și un rezultat discret, există doar un număr finit de valori posibile pentru valoarea p. Deci, desigur, nu poate avea o distribuție perfect uniformă, ceea ce implică o valoare p care ia orice valoare între 0 și 1 cu probabilitate egală. Mai jos vedem (aproximativ) câte valori posibile ia valoarea p. Cu cât dimensiunile eșantionului sunt mai mici, cu atât ne-am aștepta să vedem o divergență de la normalitate.
Și exact asta vedem. Iată distribuția completă a valorilor p pentru 100.000 de simulări ale fiecărei combinații:
… și aici este distribuția doar a valorilor p care sunt convențional „semnificative”, adică sub 0,05:
Cred că există câteva subtilități interesante aici. În special:
- Când dimensiunea celor două grupuri este o variabilă aleatorie (dar încă egală), atunci distribuția valorilor p este mult mai uniformă (dar încă nu tocmai uniformă). Practic, există multe mai multe posibilități pentru ca valorile p să le ia cu această aleatorie suplimentară.
- Valorile p sunt mai neuniforme atunci când probabilitatea de succes este 0,5 în fiecare caz, mai degrabă decât 0,2 sau 0,8.
- Valorile p pot fi încă foarte neuniforme chiar și cu o dimensiune mare a eșantionului de 1.000 de observații per grup (dacă probabilitatea de succes este aproape de 0,5)
Iată o mică perspectivă potențial interesantă asupra unuia dintre motivele pentru care acest lucru funcționează astfel – număr a diferitelor valori p unice pe care le obținem pentru fiecare combinație de dimensiune a eșantionului și probabilitățile subiacente:
size_lab `Prob=0.2` `Prob=0.5` `Prob=0.8`1 n=10 32 38 31 2 n=100 391 533 403 3 n=1000 2915 3853 2903 4 n=Pois(100) 9812 12582 9898
Nu trag nicio concluzie în nicio dezbatere de metaștiință aici. Doar notând acest fenomen interesant în distribuția valorilor p.
Iată restul codului pentru a face aceste simulări.
# Takes about 30 seconds with 100,000 reps
st <- system.time({
for(reps in c(1e3, 1e5)){
set.seed(42)
d <- expand_grid(
prob = rep(c(0.2, 0.5, 0.8), each = reps),
size = c(10, 100, 1000, NA)
) |>
mutate(size_lab = ifelse(is.na(size), "n=Pois(100)", glue("n={size}")),
prob_lab = glue("Prob={prob}")) |>
mutate(size = ifelse(is.na(size), rpois(n(), 100), size)) |>
mutate(x1 = rbinom(n(), size = size, prob = prob),
x2 = rbinom(n(), size = size, prob = prob),
# observed proportions p1 and p2 from the two populations
p1 = x1 / size,
p2 = x2 / size,
# under null hypothesis, the equal probability of both pops:
pmid = (p1 + p2) / 2,
# observed difference of the two proportions:
delt = abs(p1 - p2),
# variance of each of p1 and p2: p(1-p)/n
s1 = pmid * (1 - pmid) / size,
# standard deviation of the sum of two of those variances
sddelt = sqrt(s1 + s1)) |>
# calculate p values
mutate(pval_hand = 2 * (1 - pnorm(delt / sddelt)),
pval_fisher = tough_fisher(x1, size, x2, size),
pval_proptest = NA)
if(reps < 10000){
# when the number of reps is fairly small, I draw a pairs plot just to
# compare the different ways of calculating p values:
# - prop.test (out of the box R method). I couldn't find an easy way to vectorize this,
# is why it is only done here, via a loop. when reps is small
# - pval_hand (my home made approximate method)
# - pval_fisher (Fisher exact test, but toughened up as above to give 1 when both k1 and k2 are 0)
# The main conclusion from this is that my by-hand aproximation isn't great!
system.time({
for (i in 1:nrow(d)){
x <- d(i, )
d(i, "pval_proptest") <- prop.test(c(x$x1, x$x2), c(x$size, x$size))$p.value
}
})
# 5 seconds for reps=1000
plot1 <- function(){
d |>
select(pval_hand:pval_proptest, size_lab) |>
ggpairs() |>
print()
}
svg_png(plot1, glue("0299-pairs-{reps}"), w = 10, h = 8)
}
plot2 <- d |>
ggplot(aes(x = pval_fisher)) +
facet_grid(size_lab ~ prob_lab, scales = "free_y") +
geom_histogram(fill = "steelblue") +
scale_y_continuous(label = comma) +
labs(title = "Distribution of all p values when a null hypothesis is true",
subtitle = "Equal size binomial samples drawn from two populations with same underying probability",
x = "P value from Fisher's exact test
(when zero positive cases in both samples, p value is set to one)",
y = glue("Count of simulations (out of {comma(reps)})"))
plot3 <- d |>
filter(pval_fisher < 0.05) |>
ggplot(aes(x = pval_fisher)) +
facet_grid(size_lab ~ prob_lab, scales = "free_y") +
geom_histogram(fill = "steelblue") +
scale_y_continuous(label = comma) +
labs(title = "Distribution of significant (<0.05) p values when a null hypothesis is true",
subtitle = "Equal size binomial samples drawn from two populations with same underying probability",
x = "P value from Fisher's exact test
(when zero positive cases in both samples, p value is set to one)",
y = glue("Count of simulations (out of {comma(reps)})"))
svg_png(plot2, glue("0299-histogram-{reps}"), w = 9, h = 6)
svg_png(plot3, glue("0299-histogram-sig-only-{reps}"), w = 9, h = 6)
}
})
print(st)
# Number of unique p-values
d |>
group_by(prob_lab, size_lab) |>
summarise(number_p_values = length(unique(pval_fisher))) |>
spread(prob_lab, number_p_values)
