34

This is an elementary question, but I wasn't able to find the answer. I have two measurements: n1 events in time t1 and n2 events in time t2, both produced (say) by Poisson processes with possibly-different lambda values.

This is actually from a news article, which essentially claims that since $n_1/t_1\neq n_2/t_2$ that the two are different, but I'm not sure that the claim is valid. Suppose that the time periods were not chosen maliciously (to maximize the events in one or the other).

Can I just do a t-test, or would that not be appropriate? The number of events is too small for me to comfortably call the distributions approximately normal.

Charles
  • 1,068
  • 1
  • 7
  • 14

4 Answers4

29

To test the Poisson mean, the conditional method was proposed by Przyborowski and Wilenski (1940). The conditional distribution of X1 given X1+X2 follows a binomial distribution whose success probability is a function of the ratio two lambda. Therefore, hypothesis testing and interval estimation procedures can be readily developed from the exact methods for making inferences about the binomial success probability. There usually two methods are considered for this purpose,

  1. C-test
  2. E-test

You can find the details about these two tests in this paper. A more powerful test for comparing two Poisson means

Wazir
  • 703
  • 1
  • 6
  • 11
  • 5
    +1 Good reference, thanks. The C-test is a more rigorous version of the one I sketched, so it's well worth considering. The E-test relates a t-statistic to an appropriate distribution. Calculating that distribution involves a double infinite sum that will take $O(n_1 n_2)$ calculations to converge: fairly easy to code, probably overkill for checking the newspaper! – whuber Apr 14 '11 at 15:37
  • 2
    The author of the E-test paper wrote a simple fortran implementation to calculate p-values for two poisson means here: http://www.ucs.louisiana.edu/~kxk4695/ I ported their fortran to MATLAB here https://git.io/vNP86 – AndyL Jan 25 '18 at 21:41
12

How about:

poisson.test(c(n1, n2), c(t1, t2), alternative = c("two.sided"))

This is a test which compares the Poisson rates of 1 and 2 with each other, and gives both a p value and a 95% confidence interval.

Rob van Gemert
  • 129
  • 1
  • 2
  • 1
    It should be noted that for a two-sample problem, this uses a binomial test to compare rates – Jon Jan 07 '18 at 18:00
11

You're looking for a quick and easy check.

Under the null hypothesis that the rates (lambda values) are equal, say to $\lambda$, then you could view the two measurements as observing a single process for time $t = t_1+t_2$ and counting the events during the interval $[0, t_1]$ ($n_1$ in number) and the events during the interval $[t_1, t_1+t_2]$ ($n_2$ in number). You would estimate the rate as

$$\hat{\lambda} = \frac{n_1+n_2}{t_1+t_2}$$

and from that you can estimate the distribution of the $n_i$: they are Poisson of intensity near $t_i\hat{\lambda}$. If one or both $n_i$ are situated on tails of this distribution, most likely the claim is valid; if not, the claim may be relying on chance variation.

whuber
  • 281,159
  • 54
  • 637
  • 1,101
  • 1
    Thanks (+1), that's just the right check for this kind of off-the-cuff thing. It ended up being highly significant (p = 0.005) so the article is fine. I hope you don't mind, though, that I accepted the other answer -- it's good to know the 'real' way to do it when it matters. – Charles Apr 14 '11 at 16:08
5

I would be more interested in a confidence interval than a p value, here is a bootstrap approximation.

Calculating the lengths of the intervals first, and a check:

Lrec = as.numeric(as.Date("2010-07-01") - as.Date("2007-12-02")) # Length of recession
Lnrec = as.numeric(as.Date("2007-12-01") - as.Date("2001-12-01")) # L of non rec period
(43/Lrec)/(50/Lnrec)

[1] 2.000276

This check gives a slightly different result (100.03% increase) than the one of the publication (101% increase). Go on with the bootstrap (do it twice):

N = 100000
k=(rpois(N, 43)/Lrec)/(rpois(N, 50)/Lnrec)
c(quantile(k, c(0.025, .25, .5, .75, .975)), mean=mean(k), sd=sd(k))

     2.5%       25%       50%       75%     97.5%      mean        sd 
1.3130094 1.7338545 1.9994599 2.2871373 3.0187243 2.0415132 0.4355660 

     2.5%       25%       50%       75%     97.5%      mean        sd 
1.3130094 1.7351970 2.0013578 2.3259023 3.0173868 2.0440240 0.4349706

The 95% confidence interval of the increase is 31% to 202%.

GaBorgulya
  • 3,253
  • 15
  • 19