A fórmula para o teste Bayesiano A / B não faz sentido

10

Estou usando a fórmula do teste ab Bayesiano para calcular os resultados do teste AB usando a metodologia Bayesiana.

Pr(pB>pA)=i=0αB1B(αA+i,βB+βA)(βB+i)B(1+i,βB)B(αA,βA)

Onde

  • αA em um mais o número de sucessos de A
  • βA em um mais o número de falhas para A
  • αB em um mais o número de sucessos de B
  • βB em um mais o número de falhas para B
  • B é a função Beta

Dados de exemplo:

control: 1000 trials with 78 successes
test: 1000 trials with 100 successes

Um teste de suporte não Bayesiano padrão me fornece resultados significativos (p <10%):

prop.test(n=c(1000,1000), x=c(100,78), correct=F)

#   2-sample test for equality of proportions without continuity correction
# 
# data:  c(100, 78) out of c(1000, 1000)
# X-squared = 2.9847, df = 1, p-value = 0.08405
# alternative hypothesis: two.sided
# 95 percent confidence interval:
#  -0.0029398  0.0469398
# sample estimates:
# prop 1 prop 2 
#  0.100  0.078 

enquanto minha implementação da fórmula de Bayes (usando as explicações no link) me deu resultados muito estranhos:

# success control+1
a_control <- 78+1
# failures control+1
b_control <- 1000-78+1
# success control+1
a_test <- 100+1
# failures control+1
b_test <- 1000-100+1

is_control_better <- 0
for (i in 0:(a_test-1) ) {
  is_control_better <- is_control_better+beta(a_control+i,b_control+b_test) / 
                       (b_test+i)*beta(1+i,b_test)*beta(a_control,b_control)

}

round(is_control_better, 4)
# [1] 0

isso significa que é , o que não faz sentido, considerando esses dados.P(TEST>CONTROL)0

Alguém poderia esclarecer?

Yehoshaphat Schellekens
fonte
Uma busca de análise bayesiana com uma p-valueetiqueta? Eu pensei que os bayesianos se recusavam a ter algo a ver com valores-p.
precisa
Você está certo! apenas pensei que iria atrair mais atenção!
Yehoshaphat Schellekens 11/03/2015
@YehoshaphatSchellekens se esse foi o verdadeiro motivo de eu remover a p-valuetag, pois ela não está relacionada.
Tim
Claro, sem problemas.
Yehoshaphat Schellekens 11/03/2015

Respostas:

17

No site que você cita, há um aviso

A função beta produz números muito grandes; portanto, se você estiver obtendo valores infinitos em seu programa, trabalhe com logaritmos, como no código acima. A função log-beta da sua biblioteca padrão será útil aqui.

então sua implementação está errada. Abaixo, forneço o código corrigido:

a_A <- 78+1
b_A <- 1000-78+1
a_B <- 100+1
b_B <- 1000-100+1

total <- 0

for (i in 0:(a_B-1) ) {
  total <- total + exp(lbeta(a_A+i, b_B+b_A)
                       - log(b_B+i)
                       - lbeta(1+i, b_B)
                       - lbeta(a_A, b_A))

}

Emite total = 0,9576921, ou seja, "probabilidades de que B vencerá A no longo prazo" (citando seu link), o que parece válido, já que B no seu exemplo tem uma proporção maior. Assim, é não uma p -valor mas sim uma probabilidade de que B é maior, em seguida, A (você não esperava que fosse <0,05).

Você pode executar as simulações simples para verificar os resultados:

set.seed(123)

# does Binomial distributions with proportions
# from your data give similar estimates?

mean(rbinom(n, 1000, a_B/1000)>rbinom(n, 1000, a_A/1000))

# and does values simulated in a similar fashion to
# the model yield similar results?

fun2 <- function(n=1000) {
  pA <- rbeta(1, a_A, b_A)
  pB <- rbeta(1, a_B, b_B)
  mean(rbinom(n, 1000, pB) > rbinom(n, 1000, pA))
}

summary(replicate(1000, fun2(1000)))

Nos dois casos, a resposta é sim.


Quanto ao código, observe que o loop for é desnecessário e geralmente eles tornam as coisas mais lentas no R; portanto, você pode usar alternativamente vapplypara um código mais limpo e um pouco mais rápido:

fun <- function(i) exp(lbeta(a_A+i, b_B+b_A)
             - log(b_B+i)
             - lbeta(1+i, b_B)
             - lbeta(a_A, b_A))

sum(vapply(0:(a_B-1), fun, numeric(1)))
Tim
fonte
Hmm ... Gostaria de saber se você realmente testou a velocidade, porque vapplynão é mais vetorizado do que o forloop, pelo contrário, eles são basicamente os mesmos. Boa resposta embora.
David Arenburg 10/03
11
forLoops C / C ++ / Fortan == vetorizados; forLaço R ! = Vetorizado. Esta é basicamente a definição de vetorizado.
David Arenburg
11
@YehoshaphatSchellekens o argumento com logs não é sobre determinado software, mas sobre computação estatística geral. No exemplo no site que você cita, o código julia é fornecido - julia também é uma linguagem muito boa para programação estatística e ainda são usados ​​logs.
Tim
2
Na verdade, acabei de fazer uma pergunta sobre essa discussão exata que tivemos sobre o SO, talvez seja necessário repensar minha abordagem vapplyno futuro. Espero obter uma resposta agradável de uma vez por todas.
David Arenburg
2
Ok, depois de muito tempo pensando e resumindo os comentários e as respostas que recebi na minha pergunta, acho que tive uma compreensão geral do que vapplyrealmente é. Veja minha resposta aqui
David Arenburg