¿Cómo se calcula internamente el valor p? ¿Y cuál es el punto de datos para el que se calcula el valor p?

Question

Me he metido en una pregunta mientras repasaba la distribución de Poisson.

Aquí está:--- Una bomba nuclear falló 5 veces en 94,32 días. Dé un intervalo de confianza (IC) del 95% para la tasa de fallos por día

Así que corrí:--

poisson.test(x = 5, T = 94.32, r = 1, alternative = "two.sided", conf.level = 0.95)

Que dio una salida de:-

Exact Poisson test

data:  5 time base: 94.32
number of events = 5, time base = 94.32, p-value < 2.2e-16
alternative hypothesis: true event rate is not equal to 1
95 percent confidence interval:
 0.01721254 0.12371005
sample estimates:
event rate 
0.05301103 

Pero no he entendido en qué valor se calcula el valor p y cómo se calcula aquí?

Answer 1

Como mencionó @linog en su comentario, escribir poisson.test en R imprimirá el código fuente de la función.

La parte relevante es:

⋮
m <- r * T
PVAL <- switch(alternative, less = ppois(x, m), greater = ppois(x - 
            1, m, lower.tail = FALSE), two.sided = {
            if (m == 0) (x == 0) else {
                relErr <- 1 + 1e-07
                d <- dpois(x, r * T)
                if (x == m) 1 else if (x < m) {
                  N <- ceiling(2 * m - x)
                  while (dpois(N, m) > d) N <- 2 * N
                  i <- seq.int(from = ceiling(m), to = N)
                  y <- sum(dpois(i, m) <= d * relErr)
                  ppois(x, m) + ppois(N - y, m, lower.tail = FALSE)
                } else {
                  i <- seq.int(from = 0, to = floor(m))
                  y <- sum(dpois(i, m) <= d * relErr)
                  ppois(y - 1, m) + ppois(x - 1, m, lower.tail = FALSE)
                }
            }
        })
⋮

Y particularmente en este caso:

N <- ceiling(2 * m - x)
while (dpois(N, m) > d) N <- 2 * N
i <- seq.int(from = ceiling(m), to = N)
y <- sum(dpois(i, m) <= d * relErr)
ppois(x, m) + ppois(N - y, m, lower.tail = FALSE)

con la línea final siendo el valor p que se calcula para el x , r y T valores.

Si $X \sim \textrm{Poisson}(rT)$ el valor p para la prueba a dos bandas viene dado por: $$\sum_{x=0}^{\infty} \mathbb{P}(X = x) * \mathbb{I}_{\mathbb{P}(X=x) \le \mathbb{P}(X = 5)}$$

Donde $\mathbb{I}$ es la función indicadora.

No es más que la suma de todos los posibles valores observados que son tan probables o menos que el verdadero valor observado, dada la hipótesis nula ( $\lambda = rT$ ). Esta es la idea básica de una prueba de dos caras: queremos saber la probabilidad de observar cualquier valor tan extremo o más extremo que el que hemos observado (dado el nulo).

La razón del trabajo extra en el código fuente es porque no podemos tomar la suma sobre un número infinito de $x$ tenemos que elegir un límite después del cual las probabilidades se vuelven lo suficientemente pequeñas como para que no importen. El punto de corte depende de x , r y T .

Dirigiéndose a su comentario -

$\lambda$ (la media y la varianza de la Poisson) es igual a la tasa ( $r$ ) por la exposición ( $T$ ). Así que aquí $\lambda=94.32$ .

Para encontrar la probabilidad de dos lados, sumamos todas las probabilidades que son menores que la probabilidad del valor observado. Esto incluye las colas izquierda y derecha del Poisson. La probabilidad de la cola izquierda se puede calcular con ppois(x = 5, lambda = 94.32, lower.tail = TRUE) .

La cola derecha es un poco más complicada, ya que tenemos que encontrar el punto de la cola derecha en el que los valores se vuelven tan o más extremos (menos probables de observar) que el valor observado.

Lo siguiente funciona para estos valores, aunque necesitará un método más riguroso como el utilizado por poisson.test para trabajar con diferentes parámetros.

i <- 0:1000

# First value on the right tail that's more extreme than 5
right_tail <- min(which(dpois(i, 94.32) <= dpois(5, 94.32) & i > 5))
p_left <- ppois(5, 94.32)
# Subtract one from right_tail first to go from 0-index to 1-index of i,
# then subtract one more to include the probability of the first value (right_tail)
p_right <-  ppois(right_tail - 1 - 1, 94.32, lower.tail = FALSE)
p_val <- p_left + p_right