5
$\begingroup$

I am trying to understand the mathematics underneath the poisson.test function in the stats package in R. I am using it to compute a p-value when comparing a sample of data against another poisson rate; not another sample of data.

I've looked up the documentation around this but I cannot find anything which outlines the mathematics of the test itself, only how it is used. If anyone is able to help, or point me in the right direction I'd be greatly appreciative.

$\endgroup$
  • 1
    $\begingroup$ See the discussion here: stats.stackexchange.com/questions/350745/… which included both (one and two sample) forms of the test. It may be that you need even more detail than is there though $\endgroup$ – Glen_b -Reinstate Monica Jul 24 at 6:05
  • $\begingroup$ Questions about software are off topic here, but I think you can probably rephrase your question so that it does not involve R. $\endgroup$ – Peter Flom - Reinstate Monica Jul 24 at 11:44
  • $\begingroup$ Well I initially posted on Stack Overflow and they told me i couldn't post it there either and that I should post it here? I will refrain from posting in future. $\endgroup$ – Luke Jul 25 at 9:56
  • 4
    $\begingroup$ Some questions are not an ideal fit for any site. However, I believe your post is within our scope and I have reopened it. A question involving R (or any other statistical language) is NOT off topic if it requires statistical expertise to understand or answer (which is the case here). Indeed, R is our most popular tag by some distance; a rather bizarre occurrence if all R questions were automatically off topic. $\endgroup$ – Glen_b -Reinstate Monica Jul 28 at 1:41
7
$\begingroup$

When you say "another Poisson rate" ... if that other Poisson rate is derived from data then you are comparing data with data.

I'll assume you mean against some prespecified/theoretical rate (i.e. that you're performing a one-sample test).

You didn't state whether you were doing a one-tailed or two-tailed test. I'll discuss both

What it's doing is using the Poisson distribution with the specified rate you're testing against, and then computing the tail area "at least as extreme" (in the direction of the alternative) as the sample you got.

e.g. consider a one-tailed test; $H_0: \mu \leq 8.5$ vs $H_1: \mu > 8.5$ and the observed Poisson count of 14. Then we can compute that the upper tail at and above 14 has 0.0514 of the probability - e.g.:

Poisson upper tail probability

> 1-ppois(13,8.5)
[1] 0.05141111

(I realize this is not the best way to compute this in R - we should use the lower.tail argument instead - but wanted to make it more transparent to readers less familiar with R; by comparison ppois(13,8.5,lower.tail=FALSE) looks like an off-by-one error)

This calculation agrees with poisson.test:

> poisson.test(14,r=8.5,alt="greater")

        Exact Poisson test

data:  14 time base: 1
number of events = 14, time base = 1, p-value = 0.05141
alternative hypothesis: true event rate is greater than 8.5
95 percent confidence interval:
 8.463938      Inf
sample estimates:
event rate 
        14 

With a two-tailed test it sums those values with equal or lower probability (i.e. as with typical Fisher-style exact tests, it uses the likelihood under the null to identify what's "more extreme"):

Poisson two tailed probabilities

The probability of a 14 with Poisson mean 8.5 is about 0.024 and in the left tail the largest x-value with probability no larger occurs at 3, so the probabilities of 0,1,2 and 3 are added in:

>  1-ppois(13,8.5)+ppois(3,8.5)
[1] 0.08152019

check against the output:

> poisson.test(14,r=8.5)

        Exact Poisson test

data:  14 time base: 1
number of events = 14, time base = 1, p-value = 0.08152
alternative hypothesis: true event rate is not equal to 8.5
95 percent confidence interval:
  7.65393 23.48962
sample estimates:
event rate 
        14 

R code is publicly available -- you can check the code; in this case it bears out what I said above.

$\endgroup$
2
$\begingroup$

Glen's answer notes that you can check the code for this function, but I'm not sure if you know how to do this, so I'll augment his answer by showing you how. To check the code, just load the relevant library and type in the function name without any arguments:

library(stats)
poisson.test

function (x, T = 1, r = 1, alternative = c("two.sided", "less", 
    "greater"), conf.level = 0.95) 
{
    ...some code here...

    PVAL <- ...some code...

    ...more code here...

    structure(list(statistic = x, parameter = T, p.value = PVAL, 
        conf.int = CINT, estimate = ESTIMATE, null.value = r, 
        alternative = alternative, method = "Exact Poisson test", 
        data.name = DNAME), class = "htest")
    }
}
<bytecode: 0x0000000019efa180>
<environment: namespace:stats>

You will see from the code that the poisson.test function creates a htest object (a list that is classed as a hypothesis test) containing calculations for the test statistic, p-value, and confidence interval. The code is quite long, but a lot of it can be ignored. The parts of interest are the code to calculate the test statistic and p-value, which are about 12-15 lines of code each. You might be able to walk through it and see how each of these objects is calculated, which will tell you the mathematics they are using. This will augment Glen's answer, which confirms the output of the test in a particular case.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.