3
$\begingroup$

I'm relatively new to stats, and I have a large dataset of strings (n=114541) for which I'm interested in a binary event: whether X occurs or not (specifically, whether a specific language construct is present in a given string). X occurs rarely (443/n). X is poisson-distributed.

The strings in the dataset are arranged in a time series, and I'm trying to determine whether there's a significant change in the occurrence of X before and after a given moment in time. Because X follows a poisson distribution, I can easily calculate lambda before and after that moment in time (let's say that lambda1 = 127/35465, and lambda2 = 316/79076).

lambda1/lambda2 (or lambda1-lambda2) gives me an indication of whether the probability of X changed. What I need to know is if that change is significant (i.e., in the example above, if the difference between 0.00399 before the switch point and 0.0035 after is significant).

I've read that the C-test can test for significance between two lambda values in a poisson distribution, but I'm unsure if that's the right test. And if it is: is there a way of running it in Scipy (I've looked through the documentation, and couldn't find any mention of it)?

Improve this question
$\endgroup$
4
  • $\begingroup$ What's the basis for the assertion that "X is poisson-distributed"? $\endgroup$
    – Glen_b
    Commented Jul 26, 2014 at 2:37
  • 2
    $\begingroup$ Just a note for anyone unfamiliar with the term 'C-test' - this is just the usual 'condition on the sum and do a binomial test'. $\endgroup$
    – Glen_b
    Commented Jul 26, 2014 at 3:11
  • $\begingroup$ The nature of the data matched every description I've read of a Poisson distribution. To get some more confirmation, I plotted 114541 Bernoulli trials with p=443/114541 chance for success (which, I think, approximates my data) ... li1 = [] for i in range(0, 114541): hist(li1.append(stats.bernoulli.rvs(p=443/114541))) ... and compared this to a randomly generated Poisson distribution with lambda set to p (from above) and the size of the dataset: hist(np.random.poisson(lam=p, size=114541)) The two looked very similar, so I figured I had a Poisson distribution. Maybe wrongly so? $\endgroup$
    – zoned post meridiem
    Commented Jul 26, 2014 at 18:07
  • $\begingroup$ I'd describe that situation as "appears consistent with a Poisson distribution". $\endgroup$
    – Glen_b
    Commented Jul 26, 2014 at 22:36

2 Answers 2

Reset to default
3
$\begingroup$

At that sort of sample size and probability, the C-test should be okay.

Since this is just a binomial test, you can test it using a binomial test in Scipy. x is 127, n is 127+316 and p is 35465/(35465+79076).

There are other tests for this situation. Vanilla R offers an exact Poisson test for example.

See also:

Krishnamoorthy K. and J. Thomson (2004),
"A more powerful test for comparing two Poisson means,"
Journal of Statistical Planning and Inference, 119, pp 23–35

which indicates that an unconditional test will tend to have greater power (as is generally the case for unconditional tests in this sort of situation).

It is, of course, not exact.

In response to comments:

You're now taking $X$ to be binomial, not Poisson.

I followed your assertion that $X$ was Poisson (by which the 35465 and 79076 are simply exposure) and showed how to do the corresponding C-test.

If you want to treat the 35465 & 79076 as numbers of trials you don't need the C-test at all. You just do a straight two-sample binomial test on the trials and successes you have.

Like so (this is in R):

> prop.test(x=c(127,316),n=c(35465,79076 ))

    2-sample test for equality of proportions with continuity correction

data:  c(127, 316) out of c(35465, 79076)
X-squared = 0.9902, df = 1, p-value = 0.3197
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.0011970592  0.0003667388
sample estimates:
     prop 1      prop 2 
0.003580995 0.003996156

Incidentally, this p-value is very similar to what you get with the C-test.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I need to ask for a clarification, which is entirely due to my lack of knowledge in this area. If I want to see whether the change in # of successes after the switchpoint (316 successes, compared to 127 before, to use the above example) is significant, then wouldn't the binomial test take: x = 317, n=79076, and p=443/114541 (in other words: x = # of successes after switchpoint, n = # of trials after switchpoint, p = hypothesized probability of success across the whole time period, i.e. total successes / total trials)? $\endgroup$
    – zoned post meridiem
    Commented Jul 28, 2014 at 19:11
  • $\begingroup$ I've put my response up into my answer. $\endgroup$
    – Glen_b
    Commented Jul 28, 2014 at 23:26
  • $\begingroup$ Oops, I just noticed that I called the wrong function in R (sorry). The two sample test in R is obtained with prop.test. Fixed. (If you want a one tailed alternative, that can be done.) $\endgroup$
    – Glen_b
    Commented Jul 29, 2014 at 1:19
  • $\begingroup$ Just as an FYI for anyone else reading this: I've also found the rateratio.test package in R useful for this problem. It's designed to compare poisson ratios, and the p-values seem close to the prop.test and binom_test that Glen_b addresses above. $\endgroup$
    – zoned post meridiem
    Commented Aug 1, 2014 at 20:29
3
$\begingroup$

I wrapped Krishnamoorthy's fortran code, based on the published paper, using numpy bindings and packaged it up. Source code is on github.

Install via

pip install poisson-etest

Usage

from poisson_etest import poisson_etest

sample1_k, sample1_n = 10, 20
sample2_k, sample2_n = 15, 20
prob = poisson_etest(sample1_k, sample2_k, sample1_n, sample2_n)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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