# Simulation of maximum likelihood ratio test to test two poisson random variables

I have two random poisson variables $x_1$ and $x_2$ with value 10 and 25 respectively. I am interested to use likelihood ratio test to test the null hypothesis: $\lambda_1=\lambda_2$, versus alernate hypthesis $\lambda_1$ not equal to $\lambda_2$.

I want to use simulation to calculate power and alpha values. I would want to do it in R so any reference to R codes will be appreciated. Thanks in advance

• By "value", I guess you mean "parameter"?
– chl
Apr 10, 2011 at 10:09
– user88
Apr 11, 2011 at 14:57

This is a particularly ill-formed question.

If by "alpha" you mean Type I error, you need to go back to Square One and get definitions straight. Type I error is not something inherent in the data, or even in the hypothesis; it's a subjectively and externally applied measure of risk. And without the Type I error, you have no reference point from which to calculate Type II error, the complement of power.

Worse yet, it's not clear--after Adam's question--whether you have JUST TWO OBSERVATIONS (10 and 25), or two distributions with means of 10 and 25, and you're looking for a suitable sample size for a balanced test comparing the means. In the first case, all you can do is a likelihood ratio test that gives an approximate p-value; there's no more information to be had in two observations. In the second case, simulation can give some useful results, but you still need a value for the Type I error to get started.

• Why the new account? Apr 10, 2011 at 14:52
• You currently hold three unregistered accounts. I've merged them all onto your first account (user id=1920), but please register once and for all (use an OpenID or see the FAQ). Thanks @cardinal!
– chl
Apr 10, 2011 at 19:48
• I disagree that this question is ill formed. The 10 and 25 may be Poisson parameters, or observations, but even if they are only two observations they can be used as approximations of the parameters, so they can serve as basis of power calculation. (Alpha must be a mistake, but this does not make the question invalid.) Apr 11, 2011 at 23:54

For the simulation let's first choose sample sizes N1 and N2 for the two Poisson samples:

require(lmtest)
N1 = 20; N2 = 15


Generate a random sample and run a likelihood ratio test:

# CODE BLOCK "A"
x = c(rpois(N1, 10), rpois(N2, 25))
group = factor(c(rep('a', N1), rep('b', N2)))
m1 = glm(x ~ 1, family=poisson)
m2 = glm(x ~ group, family=poisson)
(t = lrtest(m1, m2))


Result:

Likelihood ratio test

Model 1: x ~ 1
Model 2: x ~ group
#Df  LogLik Df  Chisq Pr(>Chisq)
1   1 -158.26
2   2  -93.39  1 129.75  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Now let's run many simulations to see the power for this particular N1 and N2:

s = 1000 # 10*1000 simulations
sigs = NULL
for (i in 1:10) {
sig = 0
for (j in i:s) {
CODE BLOCK "A" COMES HERE
if (t\$Pr <= 0.05) sig = sig + 1
}
sigs = c(sigs, sig / s)
}
c(quantile(sigs, c(.025, .5, .975)), mean=mean(sigs), sd=sd(sigs))


Result:

       2.5%         50%       97.5%        mean          sd
0.991225000 0.995500000 0.999775000 0.995500000 0.003027650


Thus the power for N1 = 20; N2 = 15 is 99%.

You can calculate power for various N1 and N2 values.