# Two poisson random variables and likelihood ratio test

I have two cases:

1. Two random poisson variables $$X_1 \sim \text{Pois}(\lambda_1)$$, $$X_2 \sim \text{Pois}(\lambda_2)$$, and testing:
• Null Hypothesis: $$\lambda_1 = \lambda_2$$
• Alternate hyp: $$\lambda_1 \neq \lambda_2$$
2. Two random binomial variables $$X_1 \sim \text{Binom}(n_1, p_1)$$, $$X_2 \sim \text{Binom}(n_2, p_2)$$, where $$n_1=n_2$$:
• Null Hypothesis: $$p_1 = p_2$$
• Alternate hyp : $$p_1 \neq p_2$$

I started with likelihood ratio test for both cases, where I calculated maximum likelihood estimates for both cases, and then estimated likelihood ratio test statistics. But while reading some literature and discussion, I learned likelihood ratio test may not be a valid test in all cases, such as at low lambda values or one lambda low and other high, so for first case (poisson case), I could also explore the option of poisson random variable condition on sum of $$X_1+X_2$$ which will be $$\text{Binom}(X_1+X_2, 0.5)$$ and can be tested using Fisher exact test. Similarly for second case with two binom random variables a good choice will be McNamara test. And there are also some Bayesian tests which can be used (I do not know of any).

So I am wondering how could I do Fisher exact test for two poisson r.v condition on sum of them, and McNamara test for two binom, and if there is any bayesian test which i could use.

Preferably, i will try these test in R as I am trying to learn R.

The Bayesian test for your question is based on the integrated (rather than maximised) likelihood. So for Poisson we have:

$$\begin{array}{c|c} H_{1}:\lambda_{1}=\lambda_{2} & H_{2}:\lambda_{1}\neq\lambda_{2} \end{array}$$

Now neither hypothesis says what the parameters are, so the actual values are nuisance parameters to be integrated out with respect to their prior probabilities.

$$P(H_{1}|D,I)=P(H_{1}|I)\frac{P(D|H_{1},I)}{P(D|I)}$$

The model likelihood is given by: $$P(D|H_{1},I)=\int_{0}^{\infty} P(D,\lambda|H_{1},I)d\lambda=\int_{0}^{\infty} P(\lambda|H_{1},I)P(D|\lambda,H_{1},I)\,d\lambda$$

$$=\int_{0}^{\infty} P(\lambda|H_{1},I)\frac{\lambda^{x_1+x_2}\exp(-2\lambda)}{\Gamma(x_1+1)\Gamma(x_2+1)}\,d\lambda$$

where $P(\lambda|H_{1},I)$ is the prior for lambda. A convenient mathematical choice is the gamma prior, which gives:

$$P(D|H_{1},I)=\int_{0}^{\infty} \frac{\beta^{\alpha}}{\Gamma(\alpha)}\lambda^{\alpha-1}\exp(-\beta \lambda)\frac{\lambda^{x_1+x_2}exp(-2\lambda)}{\Gamma(x_1+1)\Gamma(x_2+1)}\,d\lambda$$ $$=\frac{\beta^{\alpha}\Gamma(x_1+x_2+\alpha)}{(2+\beta)^{x_1+x_2+\alpha}\Gamma(\alpha)\Gamma(x_1+1)\Gamma(x_2+1)}$$

And for the alternative hypothesis we have:

$$P(D|H_{2},I)=\frac{\beta_{1}^{\alpha_{1}}\beta_{2}^{\alpha_{2}}\Gamma(x_1+\alpha_{1})\Gamma(x_2+\alpha_{2})}{(1+\beta_{1})^{x_1+\alpha_{1}}(1+\beta_{2})^{x_2+\alpha_{2}}\Gamma(\alpha_{1})\Gamma(\alpha_{2})\Gamma(x_1+1)\Gamma(x_2+1)}$$

Now if we assume that all hyper-parameters are equal (not an unreasonable assumption, given that you are testing for equality), then we have an integrated likelihood ratio of:

$$\frac{P(D|H_{1},I)}{P(D|H_{2},I)}= \frac{(1+\beta)^{x_1+x_2+2\alpha}\Gamma(x_1+x_2+\alpha)\Gamma(\alpha)} {(2+\beta)^{x_1+x_2+\alpha}\beta^{\alpha}\Gamma(x_1+\alpha)\Gamma(x_2+\alpha)}$$

Which you can see that the prior information is still very important. We can't set $\alpha$ or $\beta$ equal to zero (Jeffrey's prior), or else $H_{1}$ will always be favored, regardless of the data. One way to get values for them is to specify prior estimates for $E[\lambda]$ and $E[\log(\lambda)]$ and solve for the parameters - this cannot be based on $x_1$ or $x_2$ but can be based on any other relevant information. You can also put in a few different (reasonable) values for the parameters and see what difference it makes to the conclusion. The numerical value of this statistic tells you how much the data and your prior information about the rates in each hypothesis support the hypothesis of equal rates. This explains why the likelihood ratio test is not always reliable - because it essentially ignores prior information, which is usually equivalent to specifying Jeffrey's prior. Note that you could also specify upper and lower limits for the rate parameters (this is usually not too hard to do given some common sense thinking about the real world problem). Then you would use a prior of the form:

$$p(\lambda|I)=\frac{I(L<\lambda<U)}{\log\left(\frac{U}{L}\right)\lambda}$$

And you would be left with a similar equation to that above but in terms of incomplete, instead of complete gamma functions.

For the binomial case things are much simpler, because the non-informative prior (uniform) is proper. The procedure is similar to that above, and the integrated likelihood for $H_{1}:p_{1}=p_{2}$ is given by:

$$P(D|H_{1},I)={n_1 \choose x_1}{n_2 \choose x_2}\int_{0}^{1}p^{x_1+x_2}(1-p)^{n_1+n_2-x_1-x_2}\,dp$$ $$={n_1 \choose x_1}{n_2 \choose x_2}B(x_1+x_2+1,n_1+n_2-x_1-x_2+1)$$ And similarly for $H_{2}:p_{1}\neq p_{2}$ $$P(D|H_{2},I)={n_1 \choose x_1}{n_2 \choose x_2}\int_{0}^{1}p_{1}^{x_1}p_{2}^{x_2}(1-p_{1})^{n_1-x_1}(1-p_{2})^{n_{2}-x_{2}}\,dp_{1}\,dp_{2}$$ $$={n_1 \choose x_1}{n_2 \choose x_2}B(x_1+1,n_1-x_1+1)B(x_2+1,n_2-x_2+1)$$

And so taking ratios gives:

$$\frac{P(D|H_{1},I)}{P(D|H_{2},I)}= \frac{B(x_1+x_2+1,n_1+n_2-x_1-x_2+1)} {B(x_1+1,n_1-x_1+1)B(x_2+1,n_2-x_2+1)}$$ $$=\frac{{x_1+x_2 \choose x_1}{n_1+n_2-x_1-x_2 \choose n_1-x_1}(n_1+1)(n_2+1)}{{n_1+n_2 \choose n_1}(n_1+n_2+1)}$$

And the choose functions can be calculated using the hypergeometric($r$,$n$,$R$,$N$) distribution where $N=n_1+n_2$, $R=x_1+x_2$, $n=n_1$, $r=x_1$

And this tells you how much the data support the hypothesis of equal probabilities, given that you don't have much information about which particular value this may be.

• That was a great help. Loved the derivation part. I am wondering if there is any R implementation of this test? since eventually I am going try them in R to see how they work. – Adam May 13 '11 at 23:58
• These are not complicated functions, so they do not require an "implementation", but merely evaluation of certain functions (such as hyper-geometric sampling pdf and gamma function). – probabilityislogic May 14 '11 at 4:37
• I am going to give it a try by incorporating above equations in R and see how successful I am, since I am not an expert with R yet :) – Adam May 14 '11 at 13:42

Did your reading suggest that the likelihood ratio test statistic had problems? or that the chi-squared approximation was not very good?

I expect that most of the problems are more likely the latter, the test statistic is fine, but we don't know the distribution of it under the null hypothesis. With modern computers we can estimate the distribution fairly easily and compare to that (R works great for this).

Just generate data that looks like yours but assuming the null hypothesis to be true, compute the LR statistic (or other statistic) for each simulated dataset. Now compare the same statistic for your real data to the set of simulated statistics.

Another option, depending on the nature of your data and question, is to use permutation tests where the null hypothesis is that there is no difference between your groups, they are just samples from the same population. So you mix them together and recreate the samples over and over again and compare the test statistic of interest from the original samples to the permuted samples. Fisher's exact test is actualy a permutation test.

• Yes studies shows that at some ranges of lambda parameter for poisson LRT may not be a method of choice. I have tried simulations for likelihood ratio test (for both poisson and binomial) but have not tried permutation tests yet. Based on simulation and looking at my real data, I am pretty sure likelihood ratio test will not be valid all the times. Therefore, I would like to try some other tests such as bayesian or classical (Fisher exact test McNamara) for the comparison purpose. I will also give some more thoughts to yours recommendations and will try if i could do it in R. thanks – Adam May 13 '11 at 16:44
• @Adam, Could you give a reference for the studies showing problems with the LRT ? – Scortchi - Reinstate Monica Dec 7 '12 at 1:25