As @linog mentioned in his comment, typing poisson.test
into R will print the source code for the function.
The relevant bit is:
⋮
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)
}
}
})
⋮
And particularly in this case:
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)
with the final line being the p-value that's calculated for the given x
, r
, and T
values.
If $X \sim \textrm{Poisson}(rT)$, the p-value for the two sided test is given by:
$$ \sum_{x=0}^{\infty} \mathbb{P}(X = x) * \mathbb{I}_{\mathbb{P}(X=x) \le \mathbb{P}(X = 5)}$$
Where $\mathbb{I}$ is the indicator function.
This is just the sum of all possible observed values that are as probable or less probable than the true observed value, given the null hypothesis ($\lambda = rT$). This is the basic idea of a two-sided test - we want to know the probability of observing any value as extreme or more extreme than the one we did (given the null).
The reason for the extra work in the source code is because we can't take the sum over an infinite number of $x$ values - we have to choose a cutoff after which the probabilities become small enough that they don't matter. Where this cutoff is depends on x
, r
, and T
.
Addressing your comment -
$\lambda$ (the mean & variance of the Poisson) is equal to the rate ($r$) times the exposure ($T$). So here $\lambda=94.32$.
To find the two-sided probability, we add all the probabilities that are less than the probability of the observed value. This includes both the left and right tails of the Poisson. The left tail probability can be calculated with ppois(x = 5, lambda = 94.32, lower.tail = TRUE)
.
The right tail is a bit more tricky, as we need to find the point on the right tail where the values become as or more extreme (less probable to observe) than the observed value.
The following works for these values, though you'll need a more rigorous method such as that used by poisson.test
to work with different parameters.
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